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ABSTRACT 


Numerical computations of three-dimensional time-dependent turbulent 
flows are now feasible for three main reasons. The state of the art of 
numerical methods for solving the nonlinear Navier-Stokes equations, of 
turbulent flow modeling for numerical simulations, and of fast computer 
hardware have all made three-dimensional simulation a reality. 

The present work deals with the numerical calculations of the large 
eddy structure of turbulent flows, by use of the averaged Navier-Stokes 
equations, where averages are taken over spatial regions small compared 
to the size of the computational grid. The subgrid components of motion 
are modeled by a local eddy-viscosity model. A new finite-difference 
scheme is proposed to represent the nonlinear averaged advective term 
S/SXj (U|Uj) which has fourth-order accuracy. This scheme exhibits 
several advantages over existing schemes with regard to the following: 

1. The scheme is compact— it extends only one point away in 
each direction from the point to which it is applied. 

2. It gives better resolution for high wave-number waves 
in the solution of Poisson equation. 

3. It reduces programming complexity and computation time. 

Three examples are worked out in detail. 

1. Decay of isotropic turbulence. This problem serves as a 
teri for the proposed numerical method. Comparison of 
numerical results to experimental data {Comte-Bellot and 
Corrsin, 1971) shows satisfactory agreement. This brings 
us to the conclusion that the numerical method properly 
distributes the turbulent energy among different scales, 
with proper decay rate. 

2. Homogeneous turbulent shear flow. Numerical results con- 
firm experimental data given by Champagne, et al . (1970) 
and Harris (1974) which show departure from isotropy 

and growth of length scales in the direction of shear 
as the result of the presence of mean shear. 

iii 



3. Homogeneous turbulent shear flow with system rotation. 
Numerical results show clearly the effect of rotation 
on the stability of turbulence, as was shown experi- 
mentally by Johnston (1974). 

The numerical simulation of the model equations of turbulence given in 
this work proves to be a convenient way for extracting useful information 
on the physics of a variety of turbulent flows. The numerical results 
extend existing experimental data considerably. 
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CHAPTER 1 


INTRODUCTION 


A. General Background 

In recent years, because computers have become faster and larger, 
a new approach in handling turbulent flow simulations has emerged. This 
is a three-dimensional, time-dependent numerical simulation of the large- 
scale structures of turbulent flows. In this approach, a variation of 
an old idea is used; one applies a spatial averaging operator to the 
equations of motion. The averages are taken over regions small compared 
to the size of the computational domain, and the subgrid components of 
motion are modeled with a local eddy-viscosity coefficient. 

One of the first successful attempts in using this approach was 
Dearc’orff's in a numerical study of three-dimensional turbulent channel 
flow at high Reynolds numbers (Deardorff, 1970). The main purpose of 
his work was to extend earlier work by meteorologists (Smagorinsky , 
et al . , 1965; Leith 1965), who used these methods in calculating the 
general circulation of the atmosphere in two dimensions, to fully three- 
dimensional turbulence for a laboratory problem. Other works of interest 
that followed this approach are Deardorff (1972) and Schumann (1973). 

The numerical discretization techniques that have been used in the 
mentioned works and in other investigations fall basically into two 
categories, namely, spectral methods and real -space methods. Spectral 
methods, which deal with the transformed equations in Fourier space, are 
the most accurate methods. However, when applied to the Navier-Stokes 
equations, they exhibit a major drawback. The nonlinear advective terms 
appear as convolution sums, and have to be treated in quite a complicated 
way in order to avoid aliasing errors. Another difficulty is that the 
method is capable of handling only geometrically simple problems with 
periodic boundary conditions. Detailed descriptions of the spectral 
methods can be found in the work of Orszag and co-workers {Orszag, 1969, 
1971 a,b; fox and Orszag, 1973; Orszag and Patterson, 1972). 



Real -space methods use finite differences for discretization in 
various forms. The most popular methods are second-order schemes, which 
have particular conservation and stability properties, as will be dis- 
cussed later. The main contributions to these methods can be found in the 
works of Promt (1963), Harlow and Welch (1965), Lilly (1965), Arakawa 
(1966), Williams (1969), Deardorff (1970), and others. Although finite 
difference methods are inferior to spectral methods with respect to accuracy, 
they have the advantage of not being restricted to simply shaped regions 
and periodic boundary conditions, and this generality makes them suitable 
to flows of engineering interest, i.e., flows which contain wall-bounded 
as well as boundary-free regions. 

Only a few three-dimensional time-dependent flow simulations have 
been reported so far. The main reason is the restricted capability of 
computers in the past, in terms of speed and capacity. Even today, with 
the use of the largest available computer (CDC 7600), a 32 mesh-point 
problem is the upper limit for a core-contained program, provided the 
simplest periodic boundary conditions are used, and an average simulation 
consumes about half an hour of running time. The use of peripheral 
equipment, such as tapes, discs, or other mass memory and storage, com- 
plicates the operation significantly in terms of programming and turn- 
around time. 

The reported works include those of Deardorff (1970) and Schumann 
(1973), simulating the channel flow, using 24x14x20 and 64x32x16 mesh 
points, respectively, by finite-difference methods; of Orszag and Patterson 
(1972), simulating three-dimensional homogeneous isotropic turbulence, 
using 32 mesh points by spectral methods; and of Kwak, et al. (1975), 
simulating an isotropic turbulence and a pure strained turbulent flow, 
using 32 and 16 mesh points, respectively. Clark (ongoing Ph.D. re- 
search, Stanford University) is working on turbulence modeling, using a 
64 mesh-point program. It is hoped that a 64 mesh -point problem will 
be standard for the ILLIAC IV computer, which will be operational in 
the near future. 



B. Objective of Study 

The main objective of this study is to apply new ideas, pointed 
out by Leonard (1974), to basic turbulent shear flows. According to 
this approach, the averaging operator is defined as the convolution 
integral of a normalized filtering function g(x) and the quantities 
to be averaged. This definition relates, in a clear way, the actual 
instantaneous field and the averaged field. Turbulent shear flows, which 
exist in most natural and technological flows, exhibit peculiar character- 
istics dominated by the interaction of the turbulence and the mean shear. 
When rotation is imposed on the flow, additional effects are introduced. 

It is the purpose of this study to try to understand the basic character- 
istics of such flows by means of a numerical simulation. The numerical 
method adopted in this work is the finite-difference approach in real 
space, chosen with a view of extending it in the future to problems of 
practical interest, which involve nonsimple geometry. 

In the course of study, a new fourth-order accurate finite-dif ferenc e 
scheme is proposed to represent the averaged advective term 3/3x^ (u-u^) . 
This scheme exhibits several advantages over existing schemes. It is com- 
pact, thus reducing programming complexity and computation time, and gives 
better resolution for high wave-number waves in the solution of the Poisson 
equation. In addition, it can be applied to boundaries other than 

periodic in a convenient way. The subgrid scale model is represented by 
the eddy-viscosity hypothesis with a coefficient variable in space and 
time of the form suggested by Smagorinsky (1963). The number of mesh 

g 

points used is 16 , which gives limited resolution but proved to be of 
a good size with regard to the efficient use of the CDC 7600 computer. 

Three examples are worked out in detail: 

1. Decay of isotropic turbulence. This problem serves as a 

test for the proposed numerical scheme. Comparison of numeri- 
cal results to experimental data (Comte-Bellot and Corrsin, 

1971) shows satisfactory agreement. This brings us to the 
conclusion that the numerical method properly distributes 
the turbulent energy among different wave numbers, with the 
pioper energy decay rate. 
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2. Homogeneous turbulent shear flow. Numerical results confirm 
experimental data given by Champagne, et al . (1970) and 
Harris (1974) which show departure from isotropy and growth 
of length scales in the direction of shear as the result of 
the presence of mean shear. 

3. Homogeneous turbulent shear flow with system rotation. 
Numerical results show clearly the effect of rotation on 
the stability of turbulence, as was shown experimentally 
by Johnston (1974). 

The numerical simulation of the model equations of turbulence given in 
this work proves to be a convenient way for extracting useful information 
on the physics of a variety of basic turbulent flows, and extends con- 
siderably the existing experimental data for these flows. 
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CHAPTER 2 


MATHEMATICAL FORMULATION 


The basic equations for the following analysis are the Navier-Stokes 
equations far an incompressible fluid: 

8u-(x,t) 


at 


ax, ' i V p ax. 


2 

a Uj 


V 


ax_.ax. 


i-l ,2,3 (2-1) 



( 2 - 2 ) 


where u.. , p , p , v are the velocity, pressure, density, and kinematic 
viscosity, respectively. The summation convention is implied. 

In principle, any incompressible flow is completely determined by 
the solution of these equations, provided boundary and initial conditions 
are specified. Actually, analytical solutions exist only for a few 
special cases, and one must resort to approximate solutions, by assuming 
steady-state behavior, by assuming two -dimensionality, or by dropping 
terms which are estimated to be small by an 'order of magnitude* analysis, 
etc., or one must rt^ort to numerical methods. One of the major diffi- 
culties arising in numerical calculations of turbulent flows is the wide 
range of length scales presented in the flow. Any numerical method is 
limited by its own nature, namely, it can represent only a finite range 
of length scales. The maximum wave-number that can be resolved on a 
grid of points a distance h apart is k max = ir/h, and "honest" numerical 
calculation of turbulence for high or even moderate Reynolds numbers, 
requires a number of mesh points which is far behind the capability of 
today's computers. 

As an example, consider the problem of grid-generated isotropic 

turbulence (Comte-Bellot, and Corrsin, 1971). The reported Reynolds 

number based on Taylor microscale \ and RMS turbulent velocity level 

is in the range of Re, ^60 to 70. The Kolmogorov wave-number, at which 

A -1 „ 

molecular dissipation occurs, is k^ = 34 cm . The wave-number which 
corresponds to the location of the energy peak in the three-dimensional 
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energy spectrum is - 0.5 cm" ‘ . The number of mesh points required 
to describe all scales by Fourier methods is therefore at least 

N = I[2(k K /k E )] 3 * 1.45X10 5 (2-3) 

/3 

If finite-difference methods are used, at least eight points are needed 
to describe a wave in the same order of accuracy as by the Fourier methods 

p 

(Orszag, 1969), and the number of points increases to N n» 10 . 

Our approach to overcome this difficulty is to average Eqs. (2-1) 
and (2-2). By doing so we give up the chance of resolving all scales 
of motion. Only the large eddy structure of turbulence is calculated, 
while the subgrid components of motion are modeled. 

One convenient way of averaging is space averaging over small regions 
compared to the size of the computational region. If f(x) is any function 
of position, then the averaged function f(x) is defined as the convolution 
integral of f(x) with a filtering function g(x) (Leonard, 1974): 

7(x) = gtf - /gfx-x'Jffx^dx 1 (2-4) 

*^all space 


If g is square integrable, this operator has the property 


3f _ JT 
3X.j " 3x i 


(2-5) 


(for a pro-jf see Appendix A). 

Applying definition (2-4) to Eqs. (2-1) and (2-2) and using property 
(2-5) gives ■ 






1 sp 
p a*. 


•f V 



ax. ax. 

3 3 



( 2 - 6 ) 

(2-7) 
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The difference between the local velocity and its averaged value 
if. is due mostly to the high wave-number part of the flow which cannot 
be included directly in the computations. If we set uj = u.* - ZEj and 
substitute into the term (u - u . ) , we obtain 

* V 


II. u . 

”1 Z 


u-u. + 
1 J 



( 2 - 8 ) 


The last three terms in Eq. (2-8) represent the subgrid scales and must 
be modeled in terms of the resolvable part of the velocity field, in order 
to close the equations. A convenient way of modeling these terms is by 
means of a local eddy-viscosity model. We define a stress tensor 



u ! u • + u.ul + uTuT 

i j V] i j 


(2-9) 


We further subtract the trace from a- ■ to define the stress tensor of 

I J 

the subgrid scales as 


x n = "(°ij - Ivy 1 


( 2 - 10 ) 


Substitution of Eq. (2-10) into Eq. (2-6) gives 


3U:j 

IF 


A 

3X 


. ^i^ " 3x- * 1/3 a kk^ + 3X. T ij + v 9X.3X. 

J * J J J 


3 Uj 


(2-11) 


The term 

(p7p + 1/30^) is a generalized pressure and will 

be designated 

simply as 

p* from now on. 


T ij 

is modeled by an eddy-viscosity hypothesis, 



T ij = 2Ki "ia 

(2-12) 

where 

s i j Zysx^ 3x i J 

(2-13) 


K is a proportionality constant with the dimensions ( length) 2 / time. 
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Following Smagorinsky (1963), K is defined to be proportional 
to the magnitude of the local rate-of-strain tensor of the large structure. 
To make K dimensionally consistent, a length scale should be introduced 
into the definition, and if we take it to be the mesh spacing h , then 

K - (ch) Z (2F 1a . F^) 1 ' 2 (2-14) 

Lilly (1966) estimated the model constant to be c = 0.22 by assuming 
that the wave number ir/h lies in the inertial subrange in an isotropic 
turbulence problem. However, Deardorff (1970), in his channel flow 
simulation, has found that the constant must be modified downward by 
about 505S. Later in this work, this constant will be assessed, based on 
our numerical results. 

Another model, which is similar, at least in the numerical complexity, 
to the above model is the vorticity model, where K is proportional to 
the magnitude of the local vorticity: 


where w. is given by 


K * (ch) 2 (aj 1 « i ) 1/2 

3u k 

u i " e ijk3*r 

J 


(2-15) 

(2-16) 


In a recent simulation of isotropic turbulence (Kwak, et a!., 1975), both 
models were used and no major differences were found. Based on this 
experience, we shall use Smagorinsky's model throughout this work. 

Since dissipation due to molecular interaction is much smaller than 
dissipation due to turbulence, the term 9 u-/9x-9x- appearing in Eq. 

® ij ti 

(2-6) will be excluded. To maintain incompressibility, we shall use 
a Poisson equation instead of Eq. (2-7). This equation is derived by 
taking the divergence of Eq. (2-6) and using Eq. (2-7) to eliminate some 
of the terms. 

After all these modifications the model equations are reduced to 
the following set: 
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“ 3Xj ^ u i u j^ " 9^- P * T ij “ S 
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ax.ax* 


~ ■ Qv IV * 
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(2-17) 


(2-18) 
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CHAPTER 3 


FINITE-DIFFERENCE FORMULATION 


A. General Considerations 

The major difficulty arising in the solution of the averaged set of 
equations formulated in Chapter 2 is treatment of the term S /Sx^ (u^Uj). 
One approach suggested by Leonard (1974) is first to expand u-u- in 
a Taylor series around u-tij . If a symmetrical filter is used in the 
definition 


g “7 ( Ui u j ) = ^7 f g(x-x') u.{x') UjU f ) dx* ( 3 - 1 ) 

J J j 


the result is 


(M.) - *§-(5.0.) +%-^r-^-(u,-u < ) +0(A 4 ) (3-2) 


3X. 1 j 

u 


aXj x i j / Y dx. v i j 


where A is the averaging length scale, and the constant y depends on 
the particular filtering function g(x) . For example, a Gaussian filter 
defined by 


g(x) 



exp (-6 1 x| 2 /A 2 ) 


(3-3) 


gives the value y = 24 . The second term of the RHS of Eq. (3-2) is 
the so-called "Leonard term." As we shall see, this term contributes 
substantially in controlling the proper energy transfer among the eddy 
scales, and therefore should be accounted for in the numerical calculations. 

In general, any finite-difference scheme. which is applied to the 
incompressible Navier-Stokes equations must have certain integral con- 
servation properties, if longtime integrations are required. The minimum 
requirements are mass, momentum, and energy conservation, which means that 
in the absence of dissipation and truncation errors caused by the time- 
differencing scheme, the method should have the following properties: 
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0 


(3-4) 


E(B), “i = 

p 

= 0 1= 1,2,3 (3-5) 

^ p 

If I>.u. = 0 (3-6) 

” P ’ ' 

where £ is the sum over all mesh points in an infinite (or periodic) 
domain, and (D). is the finite-difference operator for the first 
derivative in the i tn direction. Hass conservation is controlled by 
a Poisson equation, and a discussion of this property will be given in 
section C of this chapter. It can be shown easily that conditions (3-5) 
and (3-6) impose the following requirements on the advective term: 


£ (D), (u-u.) =0 i - 1,2,3 (3-7) 

P v) ' J 

£ U i (D) j (UjUj ) - 0 (3-8) 

In the past, two common second-roder schemes have been used to 
implement these requirements. The first one uses a regular mesh configura- 
tions where the variables u..,p are defined in the same location. This 
scheme is given by 


( u * u • ) - Jr (D ) . (u - u ■ ) + u * C D )-u- + u*(D Ku* 
1 it 2 i j u v o'j u j o'j u i 


0(h 2 ) 


(3-9) 


where (D ). is the central -difference operator defined by (one dimension) 

^ J 


< D o> u n = ( Vl - Vl )/2h 


(3-10) 
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The second method uses a staggered-mesh configuration with the variables 
u-» p defined as shown in Fig. 3-1. We will give the formulas for this 
case later. 

It was shown by Leonard {1974) that if the averaged equations are 
used and the term 3/ax. (u-u.) is expanded as in Eq. (3-2), the Leonard 

u u 

term is nonconservative, and preliminary estimates have shown that this 
term has a significant part in energy extraction from large scales. How- 
ever, the term 3/ax. (iLu.) is conservative. Since the Leonard is of 

2 j • j 

0(h ), a straightforward approach to retain it in a numerical calculation 
is to use a fourth-order finite-difference scheme for the term a/3Xj (u^ 
and a second-order scheme for 

2 2 

Al_3— ,-L. (5 u ) 

Y 3X £ 3X £ axj vu i u j' • 


One possible realization which maintains properties (3-7) and (3-8) for 
the term a/3Xj (iLjil..) is as follows. 

The fourth -order scheme for the term a/ ax- (u-u .) is obtained by 

J I J 

Richardson's extrapolation (Richardson, 1927). By applying the operator 
defined by Eq. (3-9) to a mesh with 2h spacing, and then combining the 
scheme derived for h with that for 2h with weights 4/3 and -1/3 , 
respectively, we have 





- 1 
3 


1 2 [^0^ j ,h ^ * u i^ D o^j,h u j + U j^ D o^j,h u i]| 

“ l|l[ (D oli,2h (u i u l i’ J + Wj,2h u j * Wj,2h u i | + ° th 


(3-11) 


2 

In order to get the Leonard term, the Laplacian a /3x £ ax £ is approximated 
by the second-order operator (D + D_).j , v/here for one dimension we have 


<»♦> u n ■ 

(“n+1 ' u n )/h 

(3-12) 

(°.) - 

^ u n ” u n-l 

(3-13) 

(D + DJ u n - 

^ u n+1 " Zu n + u n-l 

(3-14) 
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The finite-difference scheme for the Poisson equation in this realiza- 
tion is given by 

[U D o4 D o^x * U D o4 D o^y * VWVz] p = q ^ 3 “ 15 ^ 

where ^D Q is the fourth-order central -di fference operator defined by 
(one dimension) 

W u n = (' u n+2 + 8u n+l " 8u n-l * u n-2 )/12h f 3 ' 16) 

This method has been used successfully in a recent numerical simulation 
of turbulence (Kwak, et al., 1975). However, it has two major drawbacks: 

1. The scheme extends four points in the major directions 
away from the point to which it is applied. This makes 
the method inconvenient to apply in problems with boundary 
conditions that are not periodic. 

2. It gives poor resolution for high wave-numbers in the solu- 
tion of the Poisson equation, as we shall see in section 

C of this chapter. 

In the following, a new finite-difference scheme is proposed, which 
removes some of these difficulties. The scheme is compact, it uses in- 
formation from adjacent points only, and the Poisson solver follows the 
exact solution more closely. In addition, a reduction in computation 
time is achieved. 

B. The Space-Differencing Scheme 

Consider the one-dimensional advective term 

gj-tuii) = jr + t55) (3_17) 

o 9 

If we choose A /y = h /6 , we can approximate this term by the central - 
difference operator D Q with fourth-order error term: 
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(3-18) 


®r® + Tprar® ■ D o< M > + °< h4 ) 

One can interpret this result as follows. The usual central -difference 

operator, which is second-order accurate for the first derivative, is 

fourth-order accurate for the entire advective term, provided an appro- 

o 

priate filtering function is chosen. The coefficient h /6 corresponds 
to a 2h-width filter, i.e., the filter width A is twice the width of 
the mesh size h . 

This idea can be extended to higher dimensions, but the scheme turns 
out to introduce instabilities due to energy nonconservation. To overcome 
this difficulty, a staggered-mesh configuration is used. For two dimensions 
the variables are defined at the points shown in Fig. 3-1, The x-momentum 
equation is enforced on control volumes surrounding the points (i+l/2,j), 
the y-momentum equation about (i,j+l/2), and the Poisson equation for the 
pressure about the points (i,j). Using this configuration we define a 
new finite-difference operator to represent the advective term. This 
follows the idea of Eq. (3-18), i.e., the use of a lower-order finite- 
difference scheme, having a truncation error that matches the Leonard 
term, to represent the entire advective term to a higher order of accuracy. 
If f = u.u. , this operator is given by 

1 v 


8 ■?' 

8x T i,j,k 


“ a VlJ.k + TT V f 1,j+l/2,k + f i,j-l/2,k 


+ f i,j,k+l/2 + f i,j,k-l/2 ) 


(3-19) 


where 


D x f i,J,k = (f i+l/2,j,k " f i-l/2,j,k )/h 


(3-20) 


To make it clear, we expand the advective terms of the two-dimensional 
case. 
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x -momentum equation : 


ax !uu) i+1/2,j = ciD x (uu) i+l/2,j 

W m w/ 2.1 “ “V W) i*i/2.i 

.y-momentum equation ; 

&< 5s) i,j + i/2 - "'W ,1*1/2 
W (vv) i,i+1/2 “ “V vv h ,1+1/2 


+ V LD x[ (uu) i+l/2,j+l/2 + (uu) i+l/2,j-l/2] 

(3-21 ) 

+ W (K) Hl,J + < 55 >iJ (3 ’ 22) 

+ 1 ? d *Nu*i + 15 s) iJ (3 - 23) 

+ T a *V[ (VV ^+l/2.j+l/2 + <VV, i-1/2,j+l/2] 

(3-24) 


Since in the staggered-mesh configuration not all variables are available 
at points where they are needed, averages are to be taken from adjacent 
points. For example, 

^ uv) i+l/2,j+l/2 a ( u i+l/2,j * u i+l/2,j+l^ v 1,j+V2 + v i+l ,j+l/2 )/4 

(3-25) 

Using these averages in the proposed scheme, and expanding in a Taylor 
series about the points to which the scheme is applied, we obtain the 
following (two representative terms are given; other terms can be obtained 
by analogy): 
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finite-difference scheme of (uu) = 2uu^ + 

2 r 

* T? [ 6 Vxx + 4a[i xxx + 3n-«)(uu xyy 

xzz + D x D yy + °x a zz + Vxy + a z 5 xz>] + °C h4 > 

(3-26) 


+ uu 


finite-difference scheme of (uv) = uv y + u y v + 


+ £ [C3a+1 ) Cf u y v xx + | uv xxy + | u y v yy 


+ i uv + 2 u v + 3 u„ tf vj 


2 yyy yyy 


yy r 


+ I W(V« + D x 7 xy + V + 5 xx 5 y 


+ Vzz + “^zzy + 5 z D yz + S z v yz + vu zzy 


+ Vzz>] +0(h4) 


The exact expansion of the desired terms is 

£ 

y 


(3-27) 


ecu) = 2uu x * y [6U X U X X + 2K xxx + 25D xyy + 2M xzz 


+ 2u x 5 yy + 2« x u zz + 20 y u xy + 2D z u xz ] 


(3-28) 
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= uv + u v + — * u V + UV + 3u„v.,.. 

y y y L y xx xxy y yy 


+uv + u v + 3u v + 2u v + 2u v + u v 
yyy yyy yy y xy v x tu x xy u xxy v 

+ uv+uv + uv + 2v u + 2u v + vu + v u ^ I 
xx y y zz zzy z yz z yz zzy y zzj 


(3-29) 

By comparing the finite-difference scheme, Eqs. (3-26) and (3-27), to the 
exact expansion of the advective term, Eqs. (3-28) and (3-29), it can be 
seen that for a f 1 ail terms of the finite-difference scheme occur in 
the exact expansion with different multiplicative constants. In other 
words, the proposed scheme contains a Leonard term, with filtering proper- 
ties which depend on the numerical constant a . It is impossible to 
find a simple expression for this equivalent filter in terms of a , 
but one can put upper and lower bounds for A /y for any value of a . 

For example, if a - 1/3, the scheme represents a filter with 



or, in terms of the Gaussian filter, the width of the filter is between 
h and 2/2 h with most of the terms corresponding to A - /2 h . When 
a = 1, the scheme reduces to the staggered-mesh scheme first introduced 
by Fromm (1963). 

Other terms appearing in the governing equations are finite-differenced 
in the usual way, for example. 




(P« 


i+1 .3 



(3-30) 


The eddy-viscosity coefficient K is evaluated at the same place as the 
pressure p (see Fig. 3-1). 

The complete expansion of the filtered advective term according to 
Eq. (3-19) is given in Appendix D. 
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C. The Finite-Difference Poisson Equation 

The present work revealed the fact that, in the incompressible case, 
the finite-difference scheme for the Poisson equation cannot be chosen 
independently. The difference method used in the equations of motion 
dictate a particular scheme for the Poisson equation. Otherwise, a nonstable 
solution is obtained due to a rapid growth of kinetic energy. This is true 
even if a higher-order method is used for the Poisson equation. 

To exemplify this fact, consider Euler's equation. 



■ar 


i ap 

p 3X. 


(3-31) 


where d/dt is the substantial derivative. Multiply Eq. (3-31) by u^ 
and integrate over all space: 

“ - p J u 1 $7 d *. (3 ' 32 > 

Integration by parts of the RHS of Eq. (3-32) gives 



dx 

3X 1 - 



(3-33) 


Since the contribution of u.p at infinity is zero, and 3u^/3x^ - 0 
for an incompressible fluid, we have 


\ /^(Vi) dx = 0 (3-34) 

or, the kinetic energy per unit mass is conserved. The finite-difference 
form of the RHS of Eq. (3-321 is 


(3—35) 

^ P 

We use the finite-difference analog of integration by parts to expand 
this term: 


l E “iM - 




+ ?E P< D iS*) 


(3-36) 
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Each term in the RHS of Eq. (3-36) is equal to zero and energy is 
conserved. 

The use of integration by parts in a finite-difference form is con- 
sistent only if the operators of the expressions (D^p) and (D^u.) 
are identical. Therefore, a choice of an operator for the pressure- 
gradient term (D^p) in the equations of motion forces us to choose 
the same operator for the divergence term (D.u.), which brings us to 
the conclusion that the proper Poisson operator should be in the form 

a 2 

— - a n/ D-D- (3-37) 

3x!j 

To derive the Poisson operator for a staggered-mesh, one additional 
fact has to be taken into account, namely, that the variables are 
out of phase by half the mesh size relative to the pressure p , Since a 
computational cell (see Fig. 3-1) is defined as 


u(I,J,K) = u i+l/2,j,k 
v(I.J.K) = v ijj+1/2jk 

w(I,J,K) = w i f j,k+l/2 


(3-38) 


p(I,d,K) - p i( . ik 

we find that a second-order operator for 3p/3x , which appears in the 
equations of motion, is given by 


(*s\ 

\ ax /i 




[p(I+l,J,K) - p(I.J.K)] /h = (D + ) X P (3-39) 


and a second-order divergence operator for 3u/3x is given by 

(lx).. [u(I.J.K) - u(I-l,J,K)] /h = (D_) x u (3-40) 

N 7i,0,k 
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Therefore, a consistent second-order Poisson operator is given by 


(^r lx) . ^ [ P < I+1 ' J - K ' - 2pU.J,K) + p(I-1,J,K)] /h 2 = {D + DJ X P 

(3-41) 

and the entire finite-difference Poisson equation is then 

[(D + D_) x + (D + D.) y + (D + D.)J p = q (3-42) 

It is possible to derive higher-order operators for 3p/3x and 3u/3x 
compatible with each other. For example, a fourth-order Poisson operator 
is derived from the following pair: 



£-p(I+2,J,K) + 27p(I+l,J,K) 
^-u(I+l ,J,K) + 27u(I,J,K) - 


- 27p(I,J,K) + p(I-l ,K)] /24h 

(3-43) 

27u(M,J,K) + u(I-2,J,K)] /24h 

(3-44) 


and the Poisson operator is given by 


* j^p(I+3,J ,K) -54p(I+2,J,K) + 783p(I+l,J,K) - 1460p(I,J,K) 

+ 783p(I-l ,J,K) - 54p(I-2,J,K) + p(I-3,J,K)j /576h 2 (3-45) 


If we define the RHS of Eqs. (3-43) and (3-44), respectively, 
as {4 D +) x p anc * ^4 D -^ x u * ^"difference Poisson equation will 

have the symbolic form 


+ ( 4 D H D ->y + ( 4 D +4 D -) Z ] P “ 1 ' 3 ' 46) 
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In order to assess the relative accuracy of these operators, it is 

customary to compare them with the exact Fourier transform. Fourier 

transforming the various one-dimensional finite-difference Poisson operators 

2 2 

and the continuous operator 3 /3x and plotting the transformed coefficients 
versus the nondimehsional wave number kh , Fig. 3-2, one can see trie 
relative distortion introduced in the solution by the different schemes. 

The fourth-order regular-mesh curve represents Eq. (3-15). As is seen, 
the staggered-mesh configuration shows a substantial advantage over the 
regular mesh. For kh > 0.25 the difference between the fourth-order 
regular mesh and the exact Fourier method tends to grow sharply, while 
the staggered-mesh curves trace the exact solution much more closely. 

In this work, we shall use the second-order Poisson operator on 
a staggered mesh. Although inferior to its fourth-order counterpart, it 
has the advantage that the scheme extends only one point away from the 
point to which it is applied. 

The actual solution of the Poisson equation is by means of a direct 
method. The discrete Fourier transform of Eq. (3-42) is given by 


4 r.^2 mi , _.„2 mj , 2 mk"] r 

- 2 * Ism T + sin -f + sin T j p 



(3-47) 


To obtain p , we first transform q into q , then solve for p according 
to Eq. (3-47), and finally inverse transform p to obtain p . The fast 
Fourier transform algorithm (FFT) is used in the numerical realization. 

For further details of TFT methods see Cochran, et al . (1967) and 
Appendix C. 


D. The Time-Differencing Scheme 

The basic equations of motion can be reduced to the form 



s i 


(3-48) 


where s^ is a complicated nonlinear function evaluated at each time 
step for every mesh paint, by one of the methods described in the previous 
sections. 
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This equation, which is parabolic with respect to time, can be solved 
by methods which are analogous to methods for first-order ordinary 
differential equations. The existing methods can be classified according 
to the following criteria: accuracy, stability, one step or multistep, 
implicit or explicit methods. For extensive treatment of this subject 
the reader is referred to the book by Gear {1971). 

For the present problem two additional criteria were considered, 
namely, the number of evaluations of per time step, and the amount 
of storage required. These features are characteristic of problems of 
large volume and long computation time. 

Two commonly used methods, which are second-order accurate and 
require only one function evaluation per time step, are the leapfrog and 
the Adams -Bashforth methods. Both are multistep methods and require 
storage for two time steps. The leapfrog method is given by 


u.(" +1 ’ = U.'"" 15 + 2 it s>> 
1 1 1 


(3-49) 


with truncation error 


1 2 8 ^i 

- ± At — ± 

6 3t 3 


(3-50) 


and the Adams-Bashforth method is given by 
.(Ml) = u, (n) + At(3s| n > 


- s/ n_1 b/2 


(3-51) 


with truncation error 


_L At 2 ^ 
12 Tt 7 


(3-52) 


A linear-stability analysis shows that the Adams-Bashforth method is more 
stable than the leapfrog method (Lilly, 1965). 
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In this work we use the Adams -Bashforth method. The time step 
is chosen to give a Courant number of N = 0.25 where 

w 

U At 

N c = (3-53) 

J v is the magnitude of the maximum velocity appearing in the problem. 

JUtt/v 
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Fig. 3-2. Comparison of second and fourth order Poisson 
operators to the exact Fourier transform. 
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CHAPTER 4 

DECAY OF ISOTROPIC TURBULENCE 


A. General Description 

The decay of isotropic trubulence has been chosen as the first 
problem for numerical simulation. A random-velocity field which satisfies 
certain correlation and energy properties is specified within a box, and 
the evolution of this field is observed at later times. This problem, 
although not of great practical importance, has several attractive 
properties: 

1. There exist experimental data to compare with. 

2. The physical processes involved are the simplest of any 
turbulent flow. Because the mean velocity is zero pro- 
duction terms in the turbulent kinetic energy equation 
are identically zero and the only processes involved are 
energy transfer among the wave-numbers and dissipation. 
Notwithstanding, all terms in the governing equation are 
used, and can be assessed. 

3. Boundary conditions are simple. It is assumed that the 
domain of numerical integration is sufficiently large 
compared to any characteristic length associated with the 
flow field (for example, the integral scale), so that 
correlations at distances comparable to the size of the 
box are practically zero. This property allows us to 
choose convenient periodic conditions for the boundaries. 

It should be noted at this point that this problem has also inherent 
limitations. Since the experimentalists measured only certain gross 
properties, such as the shape of the three-dimensional energy spectrum 
and the rate of energy decay, it is obvious that only these properties 
can be compared with numerical results. Many of the details available 
in the numerical simulation are averaged out. It is believed, however, 
that the properties we are comparing with experimental data are the mast 
important ones for this flow, and that, if simulations predict these 
properties correctly, they represent an actual turbulent flow field 
closely enough. 
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Decay of isotropic turbulence also serves as a means of finding the 
proper value for the model constant c appearing in the subgrid scale 
model. It is chosen to fit experimental data as close as possible. 

Another method of determining c from a purely computational experiment 
is being pursued by Clark (ongoing Ph.D. research, Stanford University). 

B. The Experimental Results 

The experiment which has been chosen to be simulated is due to 
Comte-Bellot and Corrsin (1971). Measurements are given for an isotropic 
turbulent field which is generated in a wind tunnel behind a grid of 
bars. The Reynolds number based on the Taylor microscale X and RMS 
turbulent velocity level is in the range of Re^ = 60 to 70 , which 
corresponds to low levels of turbulence. Measurements are taken at three 
stations downstream in the mean flow direction. Using Taylor's hypothesis, 
the spatial loca ; ions can be regarded as temporal stations with respect 
to the evolution of homogeneous turbulence with the following relation; 

Ax = U Q At (4-1) 

where U 0 is the mean velocity. 

At these stations the three-dimensional energy spectrum E(k,t) 
is given, where E(k,t) is defined as the integral over a spherical 
shell of radius k of the trace of the spectrum tensor &.-(k) , 

* J 

which is the Fourier transform of the correlation tensor R . ^ ( r^) defined by 

^•(r) * <u 1 (x,t) Ujfx+r.t) > (4-2) 

< > is an appropriate average (see Tennekes and Lumley, 1972, Chapter 8). 
In Fourier space this function is given by 

E(k) = 2irk 2 <□. (k) u?(k)> (4-3) 

spherical shell 

It should be noted that the integral of the three-dimensional energy 
spectrum over the entire wave-number range is equal to the kinetic energy 
per unit mass 
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L 


E(k)dk = <u^ (x) u^(x) > /2 


(4-4) 


Our aim is to start the numerical simulation with initial conditions 
that correspond to the experimental data at the first station, then to 
integrate the equations of motion for the interval of time between stations 
one and two, and examine and compare numerical results with experimental 
data for the second station. Two main properties are chosen for comparison, 
namely, rate of decay of turbulent kinetic energy, and the shape of the 
three-dimensional energy spectrum. These properties are basic for this 
problem. Rate of decay shows how good our model for the subgrid scale 
dissipation is. The shape of the energy spectrum shows how accurately 
transfer and redistribution of energy among different waves are treated 
and is an indication of how well the nonlinear terms are represented by 
the numerical scheme. 

Since experimental data are regarded as representing the instantaneous 

field, while our method filters out the small-scale motion, it is required 

to filter the experimental data, so that a common basis for comparison 

is achieved. As wac shown ir Chapte 1 ' 3 the proposed numerical scheme 

2 

for the advective term *s representing a filter with two bounds on A /y . 
For a Gaussian filter we have h < A < 2/2 h with most of the terms 
corresponding to A - /2 h . If we choose the Gaussian filter with width 
of /2 h as representing the filtering process, then the three-dimensional 
energy spectrum at each station of the experiment would be 


E(k,t) = [g(k) ] 2 E(k,t) = exp(-k 2 h 2 /6)E(k,t) (4-5) 

This is the energy of the large scales o* motion. For a 16 mesh problem 
with h=2 cm, the filtered energy spectrum for the first station in the 
experiment includes about 40% of the turbulence energy. 
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C. The Initial Velocity Field 

The Fourier transform of an isotropic turbulence field of an incom- 
pressible fluid should satisfy at least four requirements: 


1. 

incompressibility 




k-jU-j (k.) 

= 0 

(4-6) 

2. 

real field 

u 1 - ( k) 

“ [ 5 i ( -y] * 

(4-7) 

3. 

isotropy 




<Uj (k) uj(k)> 

- 0 i7j 

(4-8) 

4. 

given energy spectrum 




u*(k) 

= E(k)/2Trk 2 

(4-9) 


In a numerical simulation, a modification should be made in Eq. (4-6) 
to account for the representation of derivatives by finite differences. 

In order to arrive at this modification in a simple way, consider the 
one-dimensional case. The finite-difference analog to the divergence 
operator 3/ax , which is used in the present work, is (D_) , where 

(D_) x u = (u.-ii i _ 1 )/h (4-10) 

The Fourier transforms of these operators are 

3/3x - — - > L k x (4-11) 

(D_) x -M* [(l-cos^Mt JL Sin^sij/h (4-12) 

where N is the number of mesh points in the discrete case, and m=k v /27r. 

A 

The Fourier operator associated with a continuous derivative is a pure 
imaginary number, while in the discrete case it is a complex number, which 
depends on the finite-difference scheme used. Therefore, the finite- 
difference form of the incompressibility condition should be modified to 
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k!(k) u^k) - 0 


(4-13) 


where ki is a complex vector field dependent on the finite -difference 
divergence operator D.. . This assures the numerical divergenceless 
condition in real space* 


D 1 u 1 (x) - 0 (4-14) 

A method of generating a field which satisfies Eqs. (4-7), (4-8), 

(4-9), and (4-13) is described in detail in Appendix B. The generated 
field is periodic in space, since a spectral method is used to generate it. 

D. Results and Discussion 

In Figs. 4-1 and 4-2 numerical simulation of decay of turbulent 

energy and three-dimensional energy spectra are compared to experiment 

• 3 

for the case of isotropic turbulence with mean velocity of U Q = 10 
cm/sec behind a grid of bars of diameter M = 5.08 cm. The best fit 
to the experiment was obtained with model constant c - 0.24 and 
numerical constant a - 1/3 . 

Figures 4-3 and 4-4 show the dependence of energy decay and energy 
spectrum on the numerical constant a . While energy decay differs only 
slightly for different values of a , the energy spectrum exhibits an 
irregular distribution in the a = 1 case, where energy piles up in the 
region of high wave-numbers. This case is the staggered-mesh scheme used 
by Harlow and Welch (1965), Deardorff (1970), and others. The proposed 
scheme with a = 1/3 shows a more regular behavior, which suggests that 
the Leonard term, which is included implicitly in this scheme, has an 
important effect on the correct transfer and distribution of energy, and 
should be included in the numerical scheme for proper simulation, especially 
in problems in which the energy is spread over a wide range of wave-numbers. 
This conclusion agrees with that of Kwak, et al . (1975). it is interesting 
to note that the staggered-mesh scheme with a - 1 contains implicitly 
part of the Leonard term, as can be seen from Eqs. (3-26) and (3-27), and 
was included in the past in simulations which used this scheme. 
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Sensitivity to che model constant is shown in Figs. 4-5 and 4-o. 

The effect of changing c acts mostly on the high wave-number region. 

The chosen value of c = 0.24 is 10% higher than Lilly's theoretical 
estimate of c « 0.22 (Lilly, 1966) and more than twice the value used 
by Deardorff in the channel simulation; he ran cases with c = 0.06 to 0.17 
and chose c = 0.1 to be a representative value (Deardorff, 1970). Later, 
in a paper devoted entirely to the problem of the magnitude of c , 
(Deardorff, 1971), the author suggested that the large difference between 
his numerical estimates and the value given by Lilly is due to the large- 
scale valocity shear presented in his flow. Kwak, et al. (1975) have 
found this constant to be c= .41. 

Figure 4-7 shows the skewness defined by 

S = <3 3 u/3x 3 > / <3 2 u/3x 2 > 3/2 (4-15) 

for different values of a . This function, which starts practically from 
zero, reaches a maximum value and then declines monotonically. The maximum 
value for S' is in the range based on theoretical estimates (Batchelor, 
1953, p. 172). However, T depends on the high wave-number part of the 
spectrum and our result cannot be fairly compared with experiment. 

In this chapter a numerical simulation of isotropic turbulence has 
been shown to be feasible and to give acceptable results. Computation 
time for a 16 mesh-point problem is about two minutes on the CDC 7600 
computer. In the following chapters we shall use this method with constants 
c = 0.24 and a = 1/3 to simulate homogeneous turbulence in the presence 
of shear and system rotation. 


31 



/u^XIO 



tU 0 /M-3.5 


Fig. 4-1. Rate of energy decay of isotropic turbulence, 
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FILTERED EXPERIMENT 






[q 2 / 2 ]/[ q V 2 ] 



Fig. 4-3. Dependence of rate of energy decay on 
numerical constant. 
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Fig. 4-5. Rate of energy decay for various model 
constants. 
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Fig. 4-7. Skewness for different numerical constants. 





CHAPTER 5 


HOMOGENEOUS TURBULENT SHEAR FLOW 
A. General Description 

Most natural and technological flows are shear flows. In shear 
flows most of the energy Is contained in the mean (time-averaged) velocity 
field. When the flow is turbulent, there is a coupling between the mean 
velocity and the turbulent fluctuations, and the flow characteristics are 
determined largely by this interaction. In contrast to an isotropic 
turbulent field with zero mean flow, in which turbulence levels are reduced 
by dissipation, as was shown in Chapter 4, turbulence levels in shear 
flows can be maintained or even increased by an energy -transfer mechanism 
which subtracts energy from mean flow and feeds it Into the turbulence 
in preferred directions. This effect tends to make turbulence anisotropic. 
Another characteristic of shear flows is the growth of length scales in 
the direction of shear. 

The simplest turbulent shear flow is the spatially homogeneous flow. 
Numerical results will be compared with two experiments. The first ex- 
periment is due to Champagne, Harris, and Corrsin (1970) (referred to later 
as CHC flow) and the second is due to Harris (1974). Both experiments were 
conducted in a wind tunnel. The shear flow was generated by a row of 
parallel, equal -width channels having adjustable internal resistances. 

A schematic sketch of the mean flow is shown in Fig. 5-1. In the first 
experiment, measurements were taken between the nondimensional lengths 
x/H =8.5 and x/H = 10.5 , where x is the downstream coordinate, and 
H is the height of the tunnel (H = 1 ft). The mean velocity gradient 
was r - 12.9 sec" 1 with a tunnel centerline speed of 40.7 ft/sec. 

This setup was limited in performance and left only two feet for ex- 
perimentation. In the second experiment the same setup was used, but 
the mean velocity gradient has been increased to 44 sec" 1 . If data are 
referred to a dimensionless stretched coordinate x = ( x/U ) r , the 
high- and low-shear experiments appear as a continuation of each other. 
Thus, the high-shear case can be regarded as low-shear data extended to 
distance x/H = 35. 
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The purpose of these experiments was to confirm previous hypotheses 
made by Corrsin (1963) concerning development of anisotropy with downstream 
distances, which implies asymptotic nonstationarity of the flow in a 
frame of reference convected with the mean flow. Another characteristic 
that was found in the experiments was a growth in length scales in the 
mean flow direction. The data will be presented when needed for comparison 
with calculation. 

In the numerical simulation we will try to assess these properties. 

In addition, the pressure-strain correlation tensor T . • will be computed. 
The pressure-strain term is responsible for the intercomponent energy 
transfer. According to Rotta (1951 ,1962), it destroys shear stresses, 
thus helps in redistribution of turbulence among the different components. 
Numerical results will be compared to Rotta* s estimates. 

B. Mathematical Formulation 

The mathematical representation of a spatially homogeneous turbulent 
shear flow is given by 


1 > = 
xz 

ry + u c 

(5-1) 

<V> - 

0 

(5-2) 

<W> = 

0 

(5-3) 

<P> = 

0 

(5-4) 


where U,V,W,P are instantaneous velocity components and pressure, which 
include the mean flow, r is the mean velocity gradient, < >.■ are 
averages over planes defined by ij coordinates, and < > are 
averages over the entire space. 

It is advantageous to change coordinates to a convective frame 
traveling in the x direction with centerline velocity U . The mean 

v* 

field in this frame is shown in Fig. 5-2. Let u,v,w, and p be the 
turbulence components defined as deviations from the mean flow. To 
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obtain the equations for the turbulence, we subtract the mean flow out 
by the following transformation: 



u = U ~ry 

(5-5) 


v ~ V 

(5-6) 


w » W 

(5-7) 


p = P 

(5-3) 

Substitution of these definitions in the Navier-Stokes equations 
U,V,W,P (molecular viscosity is neglected) gives 

for 

au „ 

3t 

-fH- air [(u-M-yJuj] -r & (yu) 

J 

(5-9) 

3V _ 

at 

"ply" HT (vu j ) ■ r ^T (yv) 

J 

(5-10) 

3W _ 

at 

' pH" s|r (wu j } 

J 

(5-11) 

a p 

" 3x7 axT ^ u i u j^ " 2r ax' 3x7 ( yu j) 

i J J 

(5-12) 

ax,, ax . 


In order to maintain these equations in a conservative form, no further 
simplifications should be made by means of the continuity equation. 
Using the averaging operator according to Chapter 2, we get 


H (5-13) 

J 

I - -rsf(yv) (5-H) 

J 

# - -fSf-srGSp-r&GW (5 - 15 > 

J 
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(5-16) 


1 3 2 P 


p 3x. 3X . 


3x. 3x . ^ U i U j^ “ ^ 3x 3x. ^ U j' 

i J J 


If we set u « u + u' and lump all subgrid components into the term 

3^7 T^j , we obtain 
3 


3t = ' lx " ax. ^ u+r ^) u j^ " r 3x ^ + ax. T xj t 5 ’ 17 5 

J J 


av 

at 


ly - ax'.' < v_D j> - r h + 3x7 T yj 

J J 


(5-18) 


3W 

at 


¥z “ 3x. <«J> " r 3x + 3x. T zj 

J J 


(5-19) 


2 — 

a p 

ax . 3x 1 


- 2r 9X 9x (yuj) + g 3X T 1a . 

J 1 J 


3 3 


3X. 3x . x l J 

1 u 


(5-20) 


All functionals of the form a/3x- (yu.) are treated numerically in the 

J J 

same manner as the advective term a/ax. (u.u.) , thus capturing implicitly 

J I J 

the Leonard ter.n that arises whenever a product of quantities is averaged. 
t .. is assumed to be proportional to the rate-of-strain tensor of the 

* J 

filtered turbulence components: 




id 



(5-21) 


K is given by Eq. (2-14), and we choose the model constant to be c = 0.24, 
as for the isotropic case (Chapter 4), on the assumption that the sub- 
grid scales are nearly isotropic even in the presence of mean shear. It 
should be noted that x. . now contains terms of the type uTy i.e., 

J J I 

averages containing the product of the subgrid scale turbulence and the 
mean flow. These represent interactions of the mean flow with the subgrid 
scale turbulence to produce resolvable scale fluctuations. Little is 
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known about these terms or their effects, and in using the same value of 
c as in the isotropic case we are tacitly assuming that they are 
negligible. This point deserves further investigation but will not be 
pursued here. 

a. Boundary conditions 

Since Eqs. {5-17) to (5-20) contain turbulence components only, and 
we assume that the computational domain is sufficiently large compared 
to the characteristic length scales of the flow, we shall assume u. , p 
to be periodic with respect to the computational domain. Although 
convenient for numerical implementation, this assumption has some limita- 
tions, in particular in shear flows. As we shall see, there is a growth 
of length scales in the direction of mean shear. As time goes on, these 
scales become comparable in size to the computational region, and periodicity 
assumption is questionable. In the present simulation we start with a 
field with X/H = 0.1, where A is the Taylor microscale, and H is 
the length of the computational domain. Computation is terminated when 
this ratio becomes ^0.2 . For longer simulations one must increase the 
computational box, at least in the shear direction, if periodic conditions 
are to be retained. 

b. Initial conditions 

The initial turbulence velocity field is chosen to be isotropic, 
with identical energy distribution to that of the case described in 
Chapter 4 . The turbulence level is adjusted to <u^> "^/U ^ 0.022 , 
which is approximately the initial value reported in the CHC experiment. 

c. Computational time steps 

The time steps are chosen to give a Courant number of N = 0.25 with 
respect to the mesh size h and the largest mean velocity presented 

nlaX o 

in the flow field. For U v = U/2 = 6.45 ft/sec and h = 6.25 x 10“^ ft, 

Ilia a u ^ 

we find that At - 0.25 h/U mav = 2.422 x IQ"' 3 sec. Using the Taylor 

ilia a 

hypothesis for the convective velocity U = 40.7 ft/sec, the following 
relation holds: 
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10 time steps <■ 


1 ft 


These quantities will be used throughout the numerical simulation of the 
present and the following chapter. 

C. Results and Discussion 
a. Comparison with experiments 

Figures 5-3 and 5-4 compare numerical simulation to experiments. 

The first figure compares the turbulent levels in different directions. 

As in the experiment, the v and w components, which are perpendicular 
to the shear direction, decay mono torn' cal ly with distance, while the level 
of u , the component in the shear direction, behaves differently. At 
first it decays, then it levels off and starts to grow. The growth of 
u in the CHC experiment is not clear enough, but it is shown clearly 
in the extended experiment made by Harris. Numerical results confirm 
the experimental evidence which shows that in the presence of mean shear, 
energy is extracted from the mean shear and fed into preferred turbulence 
components. This growth makes the flow nonstationary in a reference 
frame convected with the centerline velocity U . 

V 

Note that the computed results show faster growth than the experi- 
mental results. To bring the computations into agreement with experiment 
would require increasing the subgrid scale constant c which is incon- 
sistent with what Deardorff (1970) and Schumann (1973) found. Furthermore, 
the. neglected interaction between the mean flow and the subgrid scale 
fluctuations should be a production term and would tend to reduce the 
subgrid scale constant. Thus, although the computed trends agree with 
experiment, the results do differ quantitatively and the roost obvious 
possible causes would change the results in the wrong direction. The 
remaining possibility, i.e., that the difference is due to improper 
initial conditions, is also not likely to be correct as an isotropic 
initial condition should yield the smallest growth rate for a given 
turbulence level. 
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The second figure compares the Reynolds shear stress -<uv>. 

As we shall see later, this shear stress acts together with the mean 
shear to create a production term for the turbulence in the x direction. 
This term is important also in two-dimensional turbulence modeling. As 
is shown, computed results agree quite well with experiment. These com- 
parisons, although limited, bring us to the conclusion that the numerical 
simulation properly represents the main features of an actual shear flow. 

b. Additional numerical results 

In Figs. 5-5 to 5-11 additional information is extracted from the 
numerical simulation. Some of this information is either difficult or 
impossible to obtain experimentally. Examples in the first category are 
spatial correlations and length scales, and in the second category are 
pressure-strain correlations and the "return to isotropy" simulation. 

This information extends considerably the existing data set for shear 
flows and can be of help in two-dimensional turbulence modeling. A 
detailed description of these data follow. 

Figure 5-5 shows the development of total turbulent kinetic energy 
u.u 1 with time (or distance). Since the initial field is isotropic, 
energy starts to decay at an equivalent rate as in the case of isotropic 
turbulence with zero mean flow (Chapter 4). As anisotropy develops, 
enough energy is extracted from the mean fiow to compensate for dissipation, 
and to allow the net energy to increase. There is no reason to believe 
that the growth will not remain indefinitely; the steady mean field 
implies that it is being pumped in such a way as to maintain its energy. 

Figure 5-6 shows the development of turbulent levels in each direction. 
Figure 5-7 shows the three components of shear stress -<{Lu- > 

mm. U 

(normalized on the RMS values of u* and u.) . Ideally, initial conditions 

* v 

for these terms should be zero, if the field is truly isotropic. Prac- 
tically, it is difficult to achieve this requirement in a 16 problem. 
Nevertheless, the flow adjusts itself after a short period (about thirty 
steps) and the terms -<uw> and -<vw> are reduced to zero (within 
an error band of 0.01). The value of -<uv> , on the other hand, grows 
to an approximate value of 0.5. 
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Figure 5-8 shows the development of length scales A. in each 

J 

direction. The length scale in the shear direction grows considerably, 
while in the other directions there is not a significant change. Although 
experimental data for the length scales exist in the CMC experiment, 
comparison can be made qualitatively only. The reason is the large 
difference in the initial length scales of the experiment ( X .01 6 ft) 
and the numerical simulation {X^.10 ft). This difference comes mainly 
from the fact that in the numerical simulation we are using filtered 
rather than instantaneous quantities, and there is a considerable difference 
between (3u/3x) and (3u/3x) which occur in the definitions of 
A and X , respectively. 

The explanation for the observed behavior of the turbulence quantities 

and length scales that we prefer is that small nearly circular, eddies 

with their vorticity in the z direction combine to form elongated eddies 

with their major axes in the x direction. Such a process, similar to 

the vortex-pairing process put forward by the USC group (Winant and 

-2 -2 

Browand, 1974) would lead to an increase in <u > , a decrease in <v > , 

_ „ 

an increase in A v , and a constant A,. . It will also decrease <w > , 
which was observed but is not shown in the figures. 

Figure 5-9 shows nine components of spatial correlations. The most 
significant change occurs in the R^(x,0,0) component. The reason for 
that is the growth of length scales in the x direction which makes the 
flow correlated over larger distances. It should be noted that the least 
changes occure in the z direction, which means that, most of the inter- 
action takes place in the (x-y) directions. 

Figure 5-10 shows the development of the pressure-strain correlation 
tensor T^ . , where 



In an isotropic field of turbulence these terms should be identically 
zero. As can be seen, the initial values, as with the stresses - < u . u - > , 

■ J 

are not zero, and a period of adjustment should be allowed before any 
conclusion can be drawn from the simulation. In the present case we 
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shall refer to data at time step L = 100, where the T. . show regular 
behavior, and the length scales have not increased so much as to invalidate 
the periodic boundary conditions used in this simulation. More details on 
T, . will be given later in this section. 

* J 

Figure 5-11 shows a "return to isotropy" simulation from a sheared 
state. As the mean shear is turned off, the flow returns to isotropy 
in a peculiar way. When the amount of anisotropy is large, there is a 
transfer of energy between u and v which causes initially an increase 
in the v component. When the anisotropy diminishes, dissipation effects 
take over and all three components decrease toward isotropy. No such 
effect is observed in the w component, which suggests that most of the 
interaction and energy exchange is two-dimensional. This can be explined 
physically by the rotation of the elongated eddies that were formed by 
pairing. (In the absence of strain an elongated vortex will rotate.) 

c. Pressure-strain correlation 

The energy and stress balance equations for a homogeneous shear flow 
are given (Hinze, 1959, p. 252) by 


it 4 q2> 

= < -uv > r + T.ji 

+ D 

sgs 

(5-23) 

lit < ?^ 2> 

= + T 

22 

+ D 

sgs 

(5-24) 


n 

+ 

— i 
CO 

to 

+ D 

sgs 

(5-25) 

■k 

- < v > r - Tj 2 

+ D 

sgs 

(5-26) 


D$g S stands for a subgrid dissipation term which will not be specified 
explicitly here. T^ is defined in Eq, (5-22). 


Initially, when the flow is nearly isotropic, T. . is identically 

O 1 J 

zero. As time elapses, the term <vSr in Eq. (5-26), which is positive 


definite, causes a positive increase in the stress term <-uv> . This, 

-2 

in turn, increases the turbulence level of <u > in Eq. (5-23), thus making 
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the flow anisotropic. When the field is anisotropic, the terms T. • are 

I J 

no longer zero. The values of T.^ grow in a way which tends to equalize 

the velocity levels and reduce the shear stresses. Qualitatively, the 

values of T. . should be as follows: 

* J 

T .,1 c 0 , T 22 > 0 , T 33 > 0 , T-j g > 0 (5-27) 

This is indeed what we find in Fig. 5-10. 

Because of the importance of the term T. . , and the extreme difficulty 

I J 

of measuring it, it has long been a subject of model estimation. One of 
the simplest and most popular models is Rotta's estimate (1951 , 1962), 
where he assumed that the pressure-strain correlation tensor is proportional 
to the amount of anisotropy presented in the flow: 

T ij ■ - ADb ia (5 ' 28 > 

where D is the isotropic dissipation given by 

-D = ^ (q 2 /2) (5-29) 

and b.. is the anisotropy tensor defined by 

• J 

b u ■ - i z V 3)/ " z (5 - 30) 

<R ij = < 5 iV- q2 ■ R n > 

This model has been checked with data obtained from the “return to 


isotropy" 

simulation 

(Fig. 

5-11). 

The results {for L = 280) are as follows 

a. 

component 

T n 

gives a 

constant A = 2.3 

b. 

component 

T 22 

gives a 

constant A = 3.8 

c. 

component 

T 33 

gives a 

constant A = 1 .4 


These results agree with estimates made by Norris and Reynolds 
(Reynolds, 1975). They recommend a value of A = 5/2 , which is the 
average value of what we find for the separate components of T. .. 

I J 
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Fig. 5-1. Schematic sketch of the homogeneous turbulent 
shear flow. 
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Fig. 5-2. Mean velocity in a convective frame. 
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Fig. 5-3. Turbulence levels in homogeneous shear flow, 
comparison with experiment. 
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, Shear stress in homogeneous shear flow* 
comparison with experiment. 
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Fig. 5-7. Development of shear stresses in homogeneous 
shear flow. 
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Fig. 5-9. Spatial correlations in homogeneous shear flow. 






Fig. 5-10. Development of pressure-strain correlations 
in homogeneous shear flow. 
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Fi g 0 5-11. Return to isotropy from homogeneous shear. 
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CHAPTER 6 


HOMOGENEOUS TURBULENT SHEAR FLOW WITH SYSTEM ROTATION 
A. General Description 

The third problem to be simulated adds one more parameter to the 
flow field. Homogeneous turbulent shear flow is examined in the presence 
of system rotation. This problem is of interest in two important fields. 
The first one is atmospheric turbulence where the flow field is in a 
rotating frame of the earth, and the second is the internal flow of turbo- 
machines. Rotation has an effect on the stability of turbulence. It 
has been shown by experiments (Halleen and Johnston, 1967) that, depending 
on the magnitude and direction of rotation, either augmentation or sup- 
pression of turbulence intensities results. Johnston (1974) investigated 
the effects of rotation on boundary layers in turbomachine rotors and 
proposed a relation for the mixing-length ratio £/£ Q (used in two- 
dimensional turbulence modeling) which is based on the local gradient 
Richardson number defined by 

ri - . -aXLj - aa . (6-i) 

r 

in the form 

£/£ Q = 1 - 6 Ri (6-2) 

where £ is defined by 

£ = <uv> 1/2 /F (6-3) 

Based on experiments, 0 was found to be in the range of 2 to 6. 

Since our simulation represents a small region of the flow field 
away from solid boundaries or other interfaces, we do not claim that this 
simulation represents a real flow such as exhibited in an interior of a 
turbomachine; therefore, comparison with available experiments should be 
made carefully, if at all, at this stage. On the other hand, the experi- 
ments are involved and results are affected by effects other than rotation, 
as can be seen from the scattered data of the experiments (Fig. 6-6), 
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while the numerical simulation is "clean" in the sense that we are testing 
the exclusive effect of rotation on the flow. 

In the following, numerical results will be given for different 
Ri numbers and will be compared to experiment, and a simple explanation 
will be given for the stabilizing effect of rotation on the turbulence. 

B. Mathematical Formulation 

The equations of motion relative to a frame which is rotating in 
a constant angular velocity ft are given (Batchelor 1967, p. 144} by 

au , 

— + L[*VU = --VP - 2ft x U_ - ft x (ftxU) (6-4) 

V * U = 0 (6-5) 

The term -2ftxU_ is the Coriolis force, which is perpendicular to both 
IJ and ft , and -ftx(ftxU) is a centrifugal force, which can be written 
also as 


-Iv(ftxX) 2 (6-6) 

This term can be added to the pressure term -1/pVP to form a generalized 
pressure P*: 


P* = P + ~q (ftxX) 2 


In Cartesian tensor notation the equations of motion read 


au 


+ 3 (u -U * ) = -l|Ei.-2e 


at ax, 

J 


P 


i jk n j^k 


3 *1 


= 0 


Let fl - (0,0, SI), and the velocity field be defined as 


(6-7) 


( 6 - 8 ) 

(6-9) 
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<u > X2 = ry (6-10) 

<V> = 0 (6-11) 

< W > =0 (6-12) 

We define also a mean pressure 

<P> XZ = -pfiry 2 (6-13) 

To obtain the equations of turbulence, we subtract the mean values from 


the instantaneous field: 



u 

= U -ry 

(6-14) 

V 

= V 

(6-15) 

w 

= w 

(6-16) 

P 

* P* + pftr y 2 

(6-17) 


Substitution of these definitions in Navier-Stokes equations for U,V,W S P* > 
Eqs. (6-8) and (6-9), gives 


au 

at 


av 

at 


aw 

at 


a 2 p. 


ax.. ax.. 


■ p H ' axT Uu+ry)u d ] - r £ (yu) + 2 sjv (6-18) 

(6 - 19> 

- ? If - % < W V - r & (yw) < 6 - 20 > 

' 3x1" axT ^ u i u j' _ 2r F5T aj<T + 2SJ (ax ” ay) ^ 6 " 21 ' 


Following the same procedure of averaging described in Chapter 5, Eqs. 
(5-13) to (5-20), we get the conservative form of the filtered equations 
of motion relative to a rotating frame. 
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( 6 - 22 ) 


Bu 

at 


- H - w: ^ u+r ^ )u j ] + 2nv 

J Vi 


$1 

at 

3w 

at 


a 2 p 
ax. ax. 


" ax " ax, " r ax ^ * ax, T yj " 

J J 

“ at " axT ^ “ r ax' ^ + a IT T zj 

J J 

a a /- - > .p a a / - \ 4. o n /bv au \ . a 

■ 3x 1 3x, ^ u i u j^ “ r 3x 3X = ' y j ' + 2(5 ^ 3x " 3yJ 3x i 


(6-23) 

(6-24) 


3Xj T ij 


(6-25) 


Again, we use periodic conditions for the turbulent components , p . 
Initial conditions are isotropic with the same energy spectrum as for 
the case given In Chapter 5. 

C. Results 

Four different values for ft were used in the simulation: 

ft = | -1.33, -.591 , .727, 3.22j 

which correspond to gradient Richardson numbers: 

Ri = | .25, .1 , -.1 , -.25 J 

The mean velocity gradient was kept as in the case of ft = 0 in Chapter 

5, r - 12.9 Sec”^ . Figures 6-1 to 6-5 show the dependence of the total 

turbulent energy <u.u. > , the shear stress - < uv > , and the turbulence 
O 1 Jo * 1 

levels < u!j> ' on Ri , As can be seen, turbulent intensities and 
shear stresses are higher for larger rotations (or smaller Ri numbers) 
and vice versa, which means that rotation either stabilizes or destabilizes 
the turbulent field, depending on the magnitude and direction of ft . 

With regard to the turbulence levels, the most profound change occurs in 
the v component, where rotation changes its behavior entirely from a 
decaying mode for Ri - .25 to a rapidly growing mode for Ri = -.25 . 
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It is interesting to note that for Ri = -.25 , turbulence levels of u 
and v are almost identical (a simple explanation will be given later 
in this chapter) . 

Figure 6-6 shows the dependence of mixing-length ratio on the local 
gradient Richardson number, as was evaluated from experiments (Johnston, 
1974). As was mentioned before, this experiment is quite complicated 
and data are scattered, which makes any interpretation difficult, at best. 
Contrary to the experiment, which is conducted in stead-state conditions, 
our simulation exhibits nonstationarity. Nevertheless, if we analyze 
data from time step L = 100, at which time the flow has developed 
(but not too much to invalidate the periodic conditions), we find that 
3 'v 1.7 which is in the lower range of experimental results (Fig. 6-6). 

In an attempt to obtain simple correlations, we have plotted <uv> 

l/p 

and < uv > / (which is proportional to the mixing length) versus both Q 

and Ri . The best fit was obtained with <uv>^ 2 versus Ri and is 
shown in Fig. 6-7. We find 

< uv > 1/2 n, C ] - C 2 Ri (6-26) 

_9 

where ^ .5 and Cg ^ ,04t 
D. Interpretation of Results 

Simple shearing motion can be viewed as a superposition of a pure 
straining motion and a rigid rotation. When simple shear and rigid rota- 
tion are applied together, there is a value of rotation which cancels 
the rotational part of the shear and thus reduces the flow to a pure 
strain, the principal axes of which are at 45° to the (x-y) axes. The 


case of simple shearing motion in a rotating frame is similar. 

Let us 

write Eq. (6-8) in the form 



3U 

-Ter + U. tf . . 
at w j 


(6-27) 

where 

?\3x- ax. j 

J 1 * 

(6-28) 
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(6-29) 



(6-30) 


For simple shear with SU/8y - r and 
dimensions 



= (0,0, £2) we have in two 

(6-31) 
(6-32) 
(6-33) 



It is easy to see that if F/2 = 2£2 , then (J? . = -p-j and the motion 

I J I J 

reduces to pure strain with principal axes at 45°. 

The general form of the turbulent shear stress tensor for pure 
strain in the principal axes (two dimensions) is given by 


< u 2 > 


<u-u.> = 
1 3 


0 <v 2 > 


If we rotate this tensor 45°, we obtain 


1 

? 


( -2 -2 

< u > + < v > 

-2 -2 

< u > - < v > 


< u 2 > 

< u 2 > 


+ 


<v 2 >\ 

<v 2 >y 


(6—34) 


(6-35) 


with the result that the terms on the diagonal are equal. If we examine 
the case of pure shear in a rotating frame with r/2 = 2 £2 or n = 3.22 
(Ri - -.25), we find that within computational accuracy the diagonal 
elements of <u.u*> are equal, or < u > ' =<v > ' , as can be seen 

l J 

from Figs. 6-3 and 6-4. 
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This view also leads to a simple picture of turbulence production. 

If we focus strictly on the production terms, then the three most important 
equations for the turbulent velocities may be written: 


d< u > 
dt 

= -2<u 2 >^-| - 2<uv>^ 12 “ 2 < uv >R 12 

(6-36) 

d < v 2 > 
dt 

- -2 < v 2 >^22 " 2 < uv -j 2 + 2 < ijv >IR-j 2 

(6-37) 

d < uv > 
dt 

= +i n ) - (<u 2> + < v 2> )tJ ]2 

- (<v 2 >-<u 2 >)R 12 



(6-38) 


where fR^ 2 = ^12 + p 12 and we ^ ave notec * that ^12 “^21 12 = “^2 

Without loss of generality we can choose the principal axes to be along 

the coordinate directions so that 


12 ” ^ * 2^n = ~2 ^ 22 = s , 2IR12 = 1™ (6-39) 

and Eqs. (6-36) to (6-38) reduce to 

-2 

d - d “ > = s < u 2 > - r < uv > (6-40) 

-2 

_ s < v 2 > + r<uv> (6-41) 

pp =-r{ < v 2 > -<D 2 >)/ 2 (6-42) 

The kinetic energy obeys the equation 

-[“ ( < u 2 > + <v 2 >) = s(<u 2 > - <v 2 >) (6-43) 

so that there is no production of turbulent kinetic energy if the tur- 

-2 -2 

bulence is isotropic. If, however, < u > > <v > , there will be 

-2 -9 

growth caused by the fact that <u > grows faster than < v > decays. 
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Thus, if the turbulence is anisotropic, there will be production, and the 
production is in such a direction as to increase both the energy and the 
anisotropy. 

It is interesting to note the effect of rotation on production. 
Equations (6-40} to (6-42) are linear ordinary differential equations in 
< u > , < v > , and < uv> and thus have exponential solutions. Seeking 
solutions of the form e at , we find 



Thus, we see that the greatest production results when r ~ 0 , i.e., 
for pure strain (corresponding to Ri = -.25). Any rotation inhibits the 
production and, in fact, for r > s there is no production at all. 

When r = s , the growth is algebraic. 

Physically, rotation mixes the components, i.e., it transfers some 
of the u component to the v component. In so doing, it reduces the 
anisotropy and, indirectly, inhibits the production of turbulent kinetic 
energy. 
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“ <U V> / <U *■> <v 



NUMBER OF TIME STEPS 


Fig. 6-Z. Shear stresses for different gradient 
Richardson numbers. 
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Fig. 6-4. Turbulence level of v for different 
gradient Richardson numbers t 
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Fig* 6-6. Mixing length ratio 
number. Data from 








< u v > 



Fig. 6-7. <uv> 1/2 versus gradient Richardson number 

at various time steps. 
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CHAPTER 7 


SUMMARY 


A. Conclusions 

Three-dimensional time-dependent numerical simulations of simple 

turbulent flows using the filtered incompressible Navier-Stokes equations 

with the simple formulation of Smagorinsky for the subgrid scale eddy 

coefficient prove to give very reasonable results. Although the number 

3 

of degrees of freedom of 16 mesh points is limited* with respect to 
resolution, it was shown that useful information can be extracted from 
these simulations, which contain much of the physics of these flows. 
Furthermore, the approach is convenient and reasonably economical with 
regard to speed and capacity of available computers. 

The proposed numerical scheme for the filtered advective term, 
which contains a Leonard term implicitly, was shown to be important for 
proper energy transfer down the wave-number scale and is recommended for 
use in future investigations, in particular when long time integrations 
are needed and when the turbulent energy is spread over a wide range of 
wave numbers (high Reynolds numbers). 

The imposed periodic conditions for the turbulence, which have been 
used in all of our examples, have advantages as well as disadvantages. 

The advantages are twofold. First, the mathematical problem is well 
posed on the boundaries, and does not introduce errors into the domain 
of numerical integration. Second, periodic conditions allow the use of 
rapid routines in the solution of the Poisson equation. The disadvantage 
is that turbulent length scales must be kept small compared to the compu- 
tational dimension, and in problems where these scales grow, numerical 
simulation cannot be carried out for long times without invalidating 
this requirement. 

B. Recommendations 

Based on the experience gained in this work the following recommen- 
dations are proposed for future work in this field: 
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1. Improvement of initial conditions to allow better 
control on the statistics of the velocity and pressure 
fields such as stresses <u.u . > and pressure- 

* J 

strain correlations T. . . 

* J 

2. Extension of boundary conditions to allow other con- 
ditions than periodic. For example, treatment of a 
rotational -irrotational interface and a solid boundary. 

3. Extension of the computational region, at least in 
important directions such as the shear direction in 

a homogeneous shear problem, to allow longer integra- 
tion periods. This will require the use of peripheral 
equipment. 

4. Simulation of other basic flows for which experimental 
data exist such as vortex pairing and far field 
solutions of jets and wakes. 

This work is one of many efforts made in recent years the aim of 
which is to explore the basic problems of turbulent flow simulations, 
and it is believed that eventually it will bring about the solution of 
actual problems on the engineering level. 


APPENDIX A 


THE FILTERING OPERATOR AND ITS PROPERTIES 

The filtered function 7(x) is obtained as the convolution integral 
of the function f(x) with a filtering function g(x) 

7(x) = g(x) * f(x) = J g(x-x') f(x') dx' (A-l) 

all space 

x. _x' are position vectors whose components in a Cartesian system are 
x, y, z and x’ , y l , z‘ , respectively, and dx' is the volume element 
dx' dy 1 dz' . 

The only requirement on the filtering function g(x) is that it is 
Fourier transformable, which means that it should be square integrable. 

Si nee 


F 


a l k. g(k) f(k) 

(A-2) 

and 


= l k- g(k) f(k) 

(A-3) 

and also 

f [«w* a ;if] 

= l k. g(k) f(k) 

(A-4) 

then 




J-(q f) = 

3x i T ' ax i 

* f ■ 9 * SET 

(A-5) 

an important conclusion can be drawn: 



af 

3X- 

-JL f 
ax. 

(A-6) 
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The filtering function used in this work for averaging the 
Navier-Stokes equations is the Gaussian filter defined by 


g(x) 




exp (~6jx| 2 /A 2 ) 


(A-7) 


the Fourier transform of which is 


gCJl) ■ exp( - A 2 |kJ 2 /24) 


(A-8) 


Other functions can be used as well. For example, the running mean 
filter is defined by 


g(x) 


"T !*| <. box of sides A 

0 |x| > box of sides A 


its Fourier transform is given by 


g(k) 


3 

7T 

i-1 


sin (k^A/2) 
k i A/ 2 


(A-9) 


(A-10) 


This filter has been used by Deardorff (1970). For a detailed discussion 
of the properties of these filters see Kwak, et al . (1975). 
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APPENDIX B 


GENERATION OF INITIAL VELOCITY FIELD FOR ISOTROPIC 
TURBULENCE SIMULATION 


The initial velocity field for isotropic turbulence simulation 
should satisfy the following requirements (to simplify, we shall drop 


the averaging symbol ( )) 

k[(k) u.(k) = 0 (B-l ) 

u i (k) = [u i (-k}l* (8-2) 

< u^k) Uj ( k_) > - 0 i f j (B-3) 

u.{k) u*(k) = EfkVS'n-k 2 (B-4) 

where k.1 , are complex vector fields, and is the position 
vector in Fourier space. Let 

k! = (k -P r+aOc^j (B-5) 

u i = ( u i )r ^ u i ) r ( B ~ 6 ) 


Substituting the last two equations into Eqns. (B-l) and (B-4) s and 
separating real and imaginary parts, we have 


( k i )r( u i )r “ ( k iM u ih “ 0 

(B-7) 

( k i ) r( u 1 ) I " ^MWr = 0 

(B-S) 

{ u^j u j ) R + (u i -u i ) r - E(k)/27rk 2 

(B-9) 


To satisfy Eq. (B-3), the velocity components should be random. 
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The required field is generated in two stages. First, two fields 
l u i * 2 u i are 9 eneratec U wit* 1 k-j 1,631 and equal to (ki) R and ( k i ) j , 
respectively. Second, the fields are combined to give the required 
field that satisfy Eqs. {B-l ) to (B*4). 

Let us assume for the moment that we have generated these two fields; 


then 

- 0 (B-10) 

(k i ) I (2 U 1 ) = 0 (B_11) 

l u i l u i = 2 u i 2 u i s E (k)/2iTk 2 (B-l 2) 

Let the desired field be a linear combination of the two fields 

( u i )r = a l^l u i^R a 2^2 u i^R ^ B “ 13 ^ 

( u i ) I “ a 3 ^ 1 u f ^ I * a 4^2 u i^I (B-14) 

Substitute Eqs. ( B-l 3) and (B-l 4) into (B-7) to (B-9) and use Eqs. 
(B-10) and (B-ll) to get 

a 2 (k!) R ( 2 u .)R - a 3(kp i (-,u.) I = 0 (B-15) 

a T ( ^ I ( 1 u i ^ R " a 4 (k-^ ) r ( 2 u i ) I = 0 (B-16) 

a 2 ( 1 u i -j ) R + 2a 1 a 2 ( 1 u 1 2 u..) R + a|(2 u i 2 u i^R 

+ ( -j -[U.Jj + 2a 3 a 4 (^u. 2 u.j)j + 2 u ih (B-l 7) 

= E(k)/27rk 2 


80 


In order to determine the field u. , we have to determine the 
constants a-| , a 2 * a 3 , . This is done by solving the three algebraic 

equations (B-15) to (B— 1 7) in four unknowns (one is determined arbitrarily) , 
for every value of k. in the region k^ 5 0 . The other half-space 
kg > 0 is obtained by condition (B-2), which makes u. (x) real. 

The two separate fields are generated as follows: Consider a 

coordinate system in wave-number space* as shown in Fig. B-l . For any 
position vector k^ * the vector ki is specified, ki depends on the 
specific finite-difference scheme used for the divergence operator; for 
example, see Eq. (4-12). Let be a vector which lies in a plane P 
perpendicular to k! , and which makes an angle ? with a reference 
to vector t_ , which is the intersection of P and the plane formed by 
kl and e- . ? is chosen randomly in the range (0,27r) . u.u. is 

1 O n I 1 

chosen to be equal to E { k)/2Trk^ . Let n be another random number 
in the range (0,1) . Decompose u.u. into two parts such that 

( u i u 1 )r = n ^ u i u i 5 (B-18) 


(u i u i ) I = (l-n)(u.u.) 


(B-19) 


1/2 

Then decompose (u^u,. ) R ' e^ and (u-U - )j into their components 
along the vectors e_-j , e^ , e^ according to the formulae 


(u i> R ■ 
(U 2 ) R 

( U 3^R = 
< U 1>I = 

< U 2>I * 


-{cos? cos tj> cos 0 + sin 5 
(-cos? cost}) sin 0 + sin? 
(cos? sin 4 

-(cos? cos tjj cos 0 + sin? 
(-cos? cost}) sine + sin? 


sine ) )^ /2 

(B-20) 

cos 0 )(u.u i )^ 2 

(B-21 ) 


(B~22) 

sine Hu.u .) 1 / 2 

(B-23) 

cose ) (u i u i ) 

(B-24) 
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(B— 25) 


(u 3 )j - (cos? sin (j> ) Cu-u i ) 

This procedure is repeated twice for ( k ) R and (k.j)j . The sub- 
routine INICON is a realization of this method. 
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APPENDIX C 

THE DISCRETE FOURIER TRANSFORM AND ITS 
COMPUTATIONAL FORM 

The truncated sum of a Fourier expansion up to wave number N/2 
is given by 

N/2-1 

f(n) s F(m) exp(2^^ mn/N) n=0,...,N-l (C-l) 

m=-N/2 

The inverse of this transformation is given by 

, N-l 

F{m) ~ f(n) exp(-2TT^ nm/N) m=-N/2 N/2-1 (C-2) 

n n=0 

A more frequently used transform pair is given by the definition 

N-l 

f(n) = V F(m) exp(2Ti ^mn/N) n=0,...,N-l (C-3) 

m-0 

n N-l 

F(m) = i rE f(n) exp(-2TTX. nm/N) m~0,...,N-l (C-4) 

" n=0 

which is used in most fast Fourier transform (FFT) algorithms. 

Since exp is a periodic function in N , it is, in general, 
immaterial whether one set is used or the other. Both will recover 
the original function after applying a transform followed by an inverse. 
However, when discrete Fourier transforms are used to solve differential 
equations, the first pair should be used, otherwise the method will pick 
up an aliased solution which contains higher frequencies. Nevertheless, 
it is not required to give up the existing FFT codes, because there is 
a one to one correspondence between the two definitions. 

Consider that f(n), n=Q,...,N-l is given. Then, according to 
Eq. {C-2) we have 
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i N-l 

F(m) = ^{n) exp(-27r^. nm/N) m=-N/2,. . . ,N/2-l (C-5) 

n=0 


= jjr£ f(n) exp(-2irX nm'/N) 
N n=0 


m 1 


m+n for -N/2 <_ m < 0 
m for 0 < m < N/2-1 


m 1 now varies from 0 to N-l , as in Eqn. (C-4), but the ordering of 
F{m) has been changed. As an example, consider the case of N = 4 , Let 

f(n) - { f(0). f(l), f(2), f(3)} 

The discrete Fourier transform of f(n) using Eqn. (C-2) gives 

F(m) = { F(-Z), F(-1), F(0), F(l)) 

while using Eqn. (C-4) the order of F(m) is changed to 

F(m) = { F(0), F(l), F(-2), F(-l)} 

For a detailed description of the fast Fourier transform algorithms the 
reader should see Cochran, et al . (1967). 


< 


f 
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APPENDIX D 


THE COMPLETE EXPANSION OF THE FILTERED ADVECTIVE TERM 
FOR THE STAGGERED MESH CONFIGURATION 


In the following, a complete expansion of the term 3/ax,. (u^*) 
in three dimensions is given for the staggered-rnesh configuration which 
contains implicitly the Leonard term as part of its truncation error 
{see Chapter 3) (the overbar is dropped from u.. in the expansion). 

It (uu) i+l/2,j,k 


4h [ (u i+3/2,j,k + u i+l/2,j,k) ' ^ u i+l/2 s j,k T u i-l/2,j 


a 


+ u- 


j.K )2 ] 


+ Pw[ (u i+3/2,j,k * u i+3/2,j+l ,k + u i+1/2,j+1,k + u 1+1/2,j.k ) ‘ 

t 

“ ^ U i+l/2,j »k * u i+l/2,j+l ,k + u i-l/2,j+l,k + u i-l/2,j,k^ 

i 

+ (u i+3/2,j-1,k + u i+3/2,j ,k + u i+1/2,j,k + u i+l/2,j-l .k 1 ' 


‘ (u i+l/2,3-l,k + u i+l/2,j,k + u i-l/2,j,k + u i-l/2,j-l .k 5 


I 

+ (u i+3/2,j,k+l + u i+3/2,j,k + u i+l/2,j,k + “l+l^.j.k+P' 


' (u i+T/2,j,k+l + u i+l/2,j,k + u i-l/2,j,k + “l-l/i.j.k+H' 


* (u i+3/2,j,k + u i+3/2,j,k-1 + u i+l/2,j,k-l * u i+l/2,j,k ) ‘ 


" ( U 44/I 4 l/ + 




J i+l/2,j,k T u i+l/2,j,k-l + u i-l/ 2 J 3 ,k-l + u i-l/2,j,K' j ( D _t) 
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I^i+l/^.j.k 
W [ (u 1+V2,j,k 

" ( "i+1/2,j-l 

* W [ (u 1+3/2,J,k 


+ u i+l/2,j+1,k^ v i+l,;j+l/2,k + v i,j+l/2,k' 

,k + u i+l/2,j,k )(v i+l,j-1/2,k '■ v {,j-l/2,k ) ] 

+ u i+3/2,j+l ,k + u i+l/2,j+l,k + u i+l/2,j,k' v i+l,j+l/2,k 


" (u 1+3/2,j-l,k + u i+3/2,j ,k 

+ ( Vl/2,o,k + u 1+l/2,j+l .k 

' tu i+l/2,j-.1,k + u i+l/2,j,k 


+ u i+l/2,j,k + u i+V2,o-l,k )v i+l,M/2.k 
+ u i-l/2,j+l,k + u i-l/2,j,k^ v i,j+l/2,k 

+ u i-1/2.0.k + u i-l/2 > j-l,k )v T,M/2,k] 


+ fef [ (u i+l/2,j,k + Vl/2,j+1,k + u i+l/2,j,k+l + u i+l/2,j+l ,k+l 5 

• ( v i+1,j+l/2,k + v i >j+1/2,k + v i+l, 5+1/2, k+1 + v 1,j+1/2,k+l ) 

" ' u i+l/2,j-lA + u 1H/2,j,U * U i+l/2»j-l .k+1 + U 1+1/2.j,k+l^ 

* (v i+l,j-V2,k + v i,j-l/2,k + v i+1 ,j-l/2,k+l + v i ,j-l/2,k+l ) 

+ ^ u i+l/2,j,k-l + u i+l/2,j+l,k-l + u i+l/2.j>k + u i+l/2,j+l ,k^ 

’ ^ v i+l,j+l/2,k-l + v i,j+l/2,k-l + v i+1,j+1/2,k + v i,j+l/2,k' 

" lu i+l/2,a-1,k-l + u i+l/2,j,k-l + u 1+l/2,j-1 ,k + u i+1/2,5 

‘ (v 1+l,j-l/2,k-l + v i,j-l/2,k-1 + v i+l,j-l/2,k + V i .3-1/2. k^] 

(3-2) 
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^{°w)i +1 / 2 ,j,k = 

W [^ u i+1/2,j,k + u i+l/2,j,k+l^ w i+1 s j,k+l/2 + w i s j,k+l/2^ 

“ (u i+1/2,j ,k-l + u i+l/2,j s k ){w i+l,j,k-l/2 + W i,j,k~l/2 ) ] 

+ lk [ (u i+3/2,j,k+l + u H3/2,j,k + u i+l/2,j,k + u i+l/2 s j ,k+l )w i+1 


>j >k+ 1 /2 


" (u i +3/2, j ,k * u i+3/2,j ,k-l + u i+l/2,j,k-l + u 1+l/2 # J,k )w 1+1 ,J,k-l/2 
+ tu i+l/2,j,k+l + U i+l/2,j,k + U T-l/2,j,k * u i-l/2 > j J k+l Jw i,j J k+l/2 


" (u i+l/2,j,k + 


U i+l/2,j,k-l + 


u- 


i-l/2,j,k-l 


u i-l/2,j,k )w i ,j 


k-1/2 


+ 


M* [ (u i+1/2,j,k + u 1+l/2,j,k+l + u 
‘ <W 1+1,j,k+1/2 + w i,j,k+l/2 + w 

' (u i+l/2,j,k-l + u i+l/2,j,k + 11 
‘ (w 1+l,J,k-l/2 + w 1,j,k-l/Z + " 

+ (u i+1/2,j-l,k + u i+l/2,j-l ,k+1 
(w 1+l ,J-l,k+l/2 + w i, t 1-l,k+l/2 

' ( Vl/2,j-l,k-1 + “f+l/2,j-l,k 

' (w i+l,j-l, k-1/2 + w 1,j-1, k-1/2 


i+1/2,j+l,k + u i+l/2,j+l ,k+l J 
i+i ,j+l , k+1 / 2 + w i,j+l,k+l/2* 

i+1/2,j+l,k-l + ‘ i i+l/2.j+1,k ) 
1+l.J+l, k-1/2 + w i,j+l,k-l/2 ! 

+ U i+l/2,j,k + u i+l/2,j,k-M ) 

+ w 1+l,j,k+l/2 + w i,j,k+l/2 ) 

+ “l+1/2J,k-l + “1+1/2, j,k> 

+ W i+1,i, k-1/2 + w i,j,k-l/2 ) ] 


(D-3) 
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3x (vu) i,j+i/2,k 


_g 

4h 


[ (v i+l,o+l/2,k + v 1,j+l/2,k )(u 


1+1/2, J,k + u 1+1/2,j+l ,k J 


" (v i,j+l/2,k + v i-l,j+V2,k )(u i-l/2,j,k + u 1-l/2,j+l,k ) ] 

+ TOT [ (v 1+l,j+l/2,k + v i+l ,j+3/2,k + v i ,j+3/2,k + v i ,j+l/2,k )u i+l/2,j+l ,k 


~ (v. . 

» -» *1 


+ V. . 


1,J+l/2,k ‘ v i »j+3/2,k + v i~l ,j+3/2,k T v i-l ,j+l/2,k ;u i^l/2,j+l ,k 


+ V- 


X. 


+ (v. 


i + l,j-l/2,k * v i+l ,j+l/2,k T v i,j+l/2,k + v i ,j-l/2, k jU i+l/2J,k 

" (v i,j-l/2,k + v i,j+l/2,k + v i-1 ,j+1/2,k + v i‘t 1 ,j-1/2,k^ u i-l/2,j ,k] 

+ W\ [^ v i+l ,j+l/2,k + v iJ+1/2,k + v i+1 ,j+l/2,k+1 * v i,j+l/2,k+l^ 

* tu Hl/2J,k + u i+l/2,j+l,k + U i+l/2,j,k+l * u i+1/2,j+l ,R+1 ) 

' fv i,0+1/2,k + v t-l,j>V2,k + v i,j+1/2,k+l + v i-l ,j+l/2,k+l * 
^ u i~1/2,j ,k + u i-l/2,j+l ,k + u i-1/2 s a s k+l + u i -1/2, j+1 ,k+l 5 

+ ^ V i+1 »j+l/2,k-l * v i »j+l/2,k-l * v i+lJ+l/2,k * v i,j+1/2,k^ 

^ u i+l/2,j,k-l + u i+l/2,j+l ,k-l + Vl/2,j,k + u i+l/2,j+l ,k^ 

' (v i,j+l/2,k-l + v i-l ,j+l/2,k-l + v 1J+V 2,k * v i-l, 3+1/2,^ 

* < u i-l/2,j,k-l + u i-l/2,j+l,k-l * u i-l/2,j,k + u i-l/2,a+l,k ) ] 


+ v ■ . 


+ V * 


X. 


CP-4) 
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ay (vv) i,j+i/ 2 ,k 


a+3/2, k + V 1,o+l/2,k )2 * (v 1,o+l/2,l 


k * V i >j-l /2 


,S] 


+ ^f[< v 1+1 ,j. 


j+1/2,k + V i+1, 3+3/2, k + V i , 0+3/2, k + v i ,0+3/2, k J 


+ V- 


- (v, 


1+1,M/Z,k * V i+l,j+l/2,k * V i, 3+1/2, k + V 7,j-l/2,k }ii 
+ <V i *3+1/2, k * V i, 3+3/2, k + V i-1 ,j+3/2, k + V i-1 ,j+l/2,k )2 

(V l,j-l/2,k + V l, 3+1/2, k + V i-1, 3+1/2, k + V 1-l,3~l/2,k }2 


+ (v, , 


i *3+1/2, k V i,3*+3/2,k + V i, 3+1/2, k+1 + V i ,3*+3/2,k+l ^ 


(V„. i ;, + v. 4ilJO l( + v. . lm , „ + v _ j2 


1 ,J-l/ 2 ,k v i ,j+l/ 2 , k * V i ,j-l/ 2 , k +1 + v i, 


3+1/2, k+1 


+ (v, 


i »j+l/2,k-l + v i 9 j+3/2 , k— 1 + v i,j+l/2,k + v i,j+3/2,k ) 


+ v. . 


“ (v,- , 


1 *3-1/2, k-1 + v i, 3+1/2, k-1 + v i ,3-1/2, k + v i ,3+1/2, k^ ] 


(D-5) 
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« ( ™ ) i.j+l/ 2 ,k 


W[ {v i,j+l/2,k + V 1,j+l/2,k+l)^i,j,k + l/2 + w i,j + l,k + l/2> 

" (v i , 3*1/2, k-1 + V i,j+l/2,k )(w i,j,k-l/2 + w 1,j+l,k-l/2>] 

+ TSF[ (v i,j+l/2,k + V 1 ,j+3/2,k + v 1,j + 3/2,k+l + v i ( j + l/2,k + l) w T,j+1,k+l/2 

' (v i,j+l/2,k-1 + V i ,j+3/2,k-l * v i,j+3/2,k + V i ,j+l/2,k )w 1 ,j+l ,k-l/2 
+ ^ V i»j-l/2,k + v i,j+l/2,k + v i ,j+l/2,k+1 + v i,j-1/2,k+l K,j,k+l/2 
- (v i,j-l/2,k-l + v i ,j+1/2,k-1 + v i,j+l/2,k + v i,j-1/2,kH,j,k-l/2] 


+ 64 h |j v i ,j+l/ 2 ,k + v i,j+l/ 2 ,k+l + v i+l,j+l/ 2 ,k + 


v i+l ,j+ 1 / 2 ,k+l ^ 


(W 4 4 


i »j ,k+l /2 w i ,j+i ,k+l /2 4 w i+l ,j ,k+l /2 4 w i+i ,j+i ,k+l/ 2 ^ 
^ ,j+l/ 2 ,k~l 4 v i,j+l/ 2 ,k 4 v i+l ,j+l/ 2 ,k~l 4 v i+l , 0 + 1 / 2 ,^ 


(w, 


i,j,k-l /2 4 w i,j +1 ,k- 7/2 4 w i+l ,j,k-l /2 4 w i+l ,j+l ,k-l/ 2 } 


+ w_. 


+ (v. 


i-T,j+V 2 ,k 4 v i -1 ,j+l/ 2 ,k+l 4 v i,j+l/ 2 ,k 4 v i , 0 + 1/2 ,k+l ^ 


' (W 1-l,j,k+V2 + w i-1 »j+1 jk+1/2 + w 1,j,k+1/2 + w i ,J+1 ,k+l/2> 

" (v 1 - 1 ,j+l/ 2 ,k-l + v i-l,j+V 2 ,k + v 1 ,j+l/ 2 ,k-l + V 1 ,j+l/ 2 ,k> 

' (w i-l >0 .k-1/2 + w 1-U+1,k-l/2 + w i,j.k-l/2 + w i ,j+1 ,k-V2 } ] 


CD-6) 
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wu) i,j,k+l/2 = 

W [ (w i+l s j,k+l/2 + w i,j,k+l/2 )(u i+l/2,j,k+l + u i+l/2,j,k^ 

" (w i,j,k+l/2 +w i-l,j > k+l/2 )(u i-l/2,j,k+l + u i-l/2,d ,k^] 

+ TO [ (w i+l,j,k+3/2 + w i+1,j,k+1/2 * w iJ,k+l/2 + w i ,j ,k+3/2^ u i+l/2,j ,k+l 
' (w i,j,k+3/2 + w i,j,k+l/2 + w i-l ,j ,k+l/2 + w i-l ,j ,k+3/2 )L, i-l/2,j ,k+l 
* (w i+1,j,k+1/2 * w i+l »j»k-l/2 * w i,j,k-1/2 + W 1 J,k+l/2> u i+l/2,j,k 


- (w, , 


i ,j,k+l/2 W i»j,k-l/2 + w i--j ,k-l/2 * w i~l ,j,k-H/2^ u i-l/2,j,k 


64h [^ w i+l ,j ,k+1/2 w i ,j ,k+l/2 + w. +1 j +1 ,k+yz * w i } j+i s k+l/2^ 

* ( Vl/2,j,k+l * u i+l/2,j ,k + u i+l/2,j+l ,k+l + u i+1/2,j+l ,k^ 

" {w i,j,k+l/2 * w i-l ,j ,k+l/2 + w i J+l ,k+1/2 * w 1-l ,j+l ,k+l/2^ 
^ u i-l/2,j ,k+l + u i-l/2,j,k + u i-l/2,j+1,k+l + u i-1/2,j+1 ,k> 

+ (w i+1,j-l,k+l/2 * w 1,j-l,k+1/2 * w i+l ,j ,k+l/2 + w i ,j ,k+l/2^ 

* (ti i+]/2,j-1,k+l * u i+1/2,j-l,k * u i+l/2,j,k+l + u i+l/2,j,k^ 

“ (w i J j-l J k+l/2 * w i-l,j-l,k+1/2 + w i,j,k+l/2 * w i-l !,j ,k+l/2 5 
‘ ^ u i-l/2 f j-l ,k+l * u i-l/2,j-l ,k + u i-1/2,j,k+l * 
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ay^ wv ^i ,j ,k+1/2 


w[ (w i,j,k+l/2 + w 1 # j+l 1 k+l/Z^ v l,j+V2 l k * v i »j+l/2,k+li- 
“ (w i,j-l,k+1/2 + w i,j,k+l/2^ v i,j-l/2,k + V 1 * j-l/2,fc+l 
TW [ (w i\j,k+l/2 + w i,j‘+l,k+l/2 + w i „j+l ,k+3/2 + w i ,j ,k+3/2 )v i ,j+1/2, 


k+1 


" ^ W i »j-l »k+l/2 + W 1 ,j,k+l/2 + w i jjjk+3/2 + w i ,j-1 ,k+3/2 Jv i ,j-1/2,k+l 
+ ^ W i jk-1/2 + w i ,j+l ,k-l/2 + w i ,j+1 ,k+l/2 * w i ,j,k+l/2^ v i ,j+l/2,k 
" fw i,j-l,k-1/2 + w i,j,k-1/2 +w i # J,k+l/2 * w i ,j-l ,k+l/2 )v i ,j-l/2,k] 
*" 64h [ (w i,j,k+l/2 + W 1 ,j+l ,k+l/2 + w i+U,k+l/2 * w i+l ,j+l ,k+l/2 5 

^ V i,d+l/2,k * v i,j+1/2,k+l + Vl,j+t/2,k * v i+l ,j+l/2,k+l' 

" (w i,j-1,k+V2 + W i ,j,k+1/2 + w i+l ,j-l ,k+l/2 * w i+l,j,k+l/2 } 

' (v i,j-1/2,k + v i ,j-l/2,k+'I * v i+1J-l/2,k + v i+l ,j-l/2,k+l J 

+ (w i-l,j,k+l/2 + w i-1,j+l,k+l/2 + w i,j,k+l/2 + w i J+l ,k+1/2^ 

’ Cv i“l»0+V2,k + v j-l »j+l/2,k+1 + v i,j+1/2,k + V 1 ,j+l/2,k+l ) 

" Cw i-l,j-l,k+l/2 + w i-l ,j ,k+l/2 + w i,j-l ,k+l/2 + w i,j,k+1/2 } 

* tv i-l.J-l/2,k + v i-l ,d-l/2,k+l + v i,j-l/2,k + v i ,j-l/2,k+l 5 ] 

(D— 8) 
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4 ° h L <W i .J.k+l/2 + W i,j,k+3/2 )2 ' < w 1,j,k-l/2 - w i,j,k + 1/2 


> 2 ] 


1 -a 


* M [ (w i,j,k+l/2 + w i , j,k+3/2 + w 1+l,j,k+l/2 + W i+l,j,k + 3/2 )2 


^ i ,j ,k-]/2 w i,j,k+l/2 + w i+1 ,j,M/2 * w t+l ,j .k+1/2^ 


^ i"l*j,k+l/2 w i-1 ,j,k+3/2 + w i,j,k+l/2 + w i,j,k+3/2^ 2 
i-l>j,k-l/2 w i-l ,j,k+l/2 + w i,j,k-l/2 + w i ,j,k+l/2^ Z 
T (W 1 J.k-H/2 + w i ,j ,k+3/2 + w i ,j+l ,k+l/2 * w i , j+1 ,k+3/2 )2 
“ (w t.J,k-l/2 + w 1,J,k+l/2 + w i ,j+l ,k-1/2 + w i ,j+l ,k+T/2 )2 
^ 'i»j-l,k+l/2 + w i *j-l 3 k+3/2 + w f,j,k+l/2 + w i,j,k+3/2^ 2 


“ fa* 4 


1»j-l*k-l/2 + w i »j~l jk+1/2 + w i a j , k-1 /2 + w i,j%k+l/2^ j 


(D-9) 
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APPENDIX E 


PROGRAM TFC { I N PUT, OUTPUT, TAPE5 = I NPUT,TAPE6=OUTPUT,TAPE8,TAPE9) 
******************************************************************* 
* LANGUAGE: FORTRAN COMPUTER : CDC7G0Q COMPI LER:RUN76, LBL BERKELEY* 


* NOMENCLATURE: * 

* H MESH SIZE * 

* DT TIME STEP * 

* UO MAXIMUM MEAN VELOCITY IN SHEAR FLOW * 

* OMEGA ROTATION * 

* C TURDULENT MODEL CONSTANT * 

* FDC FINITE DIFFERENCE CONSTANT * 

* L NUMBER OF TIME STEPS * 


LARGE PR1(16,16, 16 ) , PI 1 ( 16, 1G, 16) /WAVE ( 1C, 16, 16) 

LARGE Ul( 16, 16,19), V1(16,1G, 19 ) , U1 ( 16, 16,19) 

LARGE RXK16, 16, 16 ), RYK16 , 16, 1G),RZ1C1G, 16,16) 
DIMENSION UC 22, 22, 7), VC 22, 22, 7),WC22, 22,7) 

DIMENSION EVC 20, 20,5), PRC 20, 20, 5) 

EQU I VALENC E C EVC 1, 1, 1 ) , PRC 1, 1, 1) ) 

DIMENSION SXC1S,1S,3),SY(18,18,3),SZ(18,13,3) 

DIMENSION QRC16,1G), at (16,16) 

DIMENSION NC3),M1C6),M2CS),TRRC1G,3),TRI C 16, 3), EC IS) 

CO MMO N / S CM / Q R , Q I , N,M1,M2, TRR, TR I , E, ! SI GN,PA! 

REAL LAMBDA 
PA 1 =3.141532653589793 
PA1 2=PA1*PAI 
H=l./16, 

H2»H*H 
DT=0 . 025 
C=0 . 24 
U0=0. 0 
OMEGA=0.0 
FDC=l./3. 

C1=FDC/H 

C2-C l."FDC)/4./H 
C3=C*C*H 
C4=Cl/4. 

C5=C2/1G. 

C6=C2/4, 

C7=l . / CH2*4. ) 

C8=1./C2.*H) 

C3=l./H 
Cll=Cl/2. 

C12«C2/'4. 

CALL TRRTRI 
-INITIAL DATA: 

-SUUROUT I IJE INI COM GENERATES THE INITIAL VELOCITY FIELD 
-ACCORDING TO APPENDIX B 
CALL INICON 
DO 11 K-l, 1G 
DO 11 J=l, 16 
DO 11 1=1,16 
P.X1C I , *J,K) =0, Q 
RY1C I ,0,10-0. 0 
P.Z1C l , 0, K) =0 , Q 
11 CONTINUE 

C -FIRST STEP IS ADVANCED BY EULER METHOD (ALPilA=1.0), 
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C ADAMS “BASHFQRTH METHOD IS USER THEREAFTER CALPHA-1.5) 

ALPHA-1- 0 
L=0 

usqr=o.o 

DUDX2=0. 0 
DUDX5-0. 0 
ENRG=0.0 
Dl.V-Q.O 
6 M“1 

DO 2 K=l,3 

KK=K+1B 

DO 2 J=l,16 

DO 2 1=1,16 

U1CI, J,KK)=U1U, J,K) 

V1CI,J,KK)=V1CI,J,K) 

Win, J,UiO=Wl(l, d,K) 

2 CONTINUE 

WRI TEC G, 102) L 

102 FORMATClHl, 10X, 22HNUMBER OF TIME STEP L=, 13) 

WRI TEC 6, 107) 

107 FORMATC1II ,10x, * U COMPONENT OF VEL FIELD *) 

SMALL I NCQRC 1,1), IJ1( 1,1, 9), 256) 

DO 44 1=1,16,8 
WR 1 TEC G, 106) I 

WRI TEC 6, 105) CQR C I , J) , J=l, 16) 

44 CONTINUE 

WR t TEC 6, 108 ) 

108 FORMATC 1H ,1QX,* V COMPONENT OF VEL FI ELI) *} 

SMALL I NCQRC 1, 1) , VIC 1, 1, 9), 256) 

DO 45 1=1,16,2 
WR I T E (- 6 , 1 0 6 ) I 

WR I TEC 6, 105 ) (QRC I , 0) , J=l, 16 ) 

45 CONTINUE 
WR1TEC6, 109} 

109 FORMATC lil ,10X,* W COMPONENT OF VEL FIELD *) 

SMALL I NCQRC 1, 1) , W1C 1, 1, 9), 256) 

DO 46 1=1,16,8 
WR! TEC 6, 106) I 

WRI TEC 6, 105) (QRCI, J),J»1,1G) 

46 CONTINUE 

1 IFCM.EO-1) GO TO 8 
K»7 

GO TO 16 

8 DO 20 K5=l, 7 

K=K5 

K2=K5+13 

IFCK2.GT.16) K2=K2-16 
16 SMALL I NCQRC 1, 1),U1C1, 1, K2) , 256) 

DO 29 J=l, 16 
DO 29 1=1,16 
JJ=J+3 
11=1+3 

un i , jj, io =qrc i , j) 

29 CONTINUE 

SMALLI NCQRC 1,1), VIC 1, 1,K2), 256) 


OKIffMAL PAGK !g 
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DO 83 1=1,16 
DO 83 J=1,1G 
JJ-J+3 
I 1=1 + 3 

V(II,JJ,K) =QR( 1,0) 

89 CONTINUE 

5MALLI N(QR(1, 1),W1(1, 1,1C2), 25G) 

DO 99 J=l,16 
00 99 1=1,16 
Jij=J + 3 
] I =1+3 

I/O I, JJ,tO*QR(|,d) 

99 CONTINUE 

c EXTEND VELOCITY FIELD PERIODICALLY IN X, Y DIRECTIONS 

DO 24 J=4, 13 
UC1, J,I0=U(17, J,K) 

V{1, J,K)=V(17, J,K) 
i/Cl, J,K)=W(17, J,K) 

U(2, J,iO=U(18, J,K) 

VC 2, J, 1C) =Vt 18, J, K) 

WC2, J,K)=W£1S; J,K) 

U{ 3, <J,iO =U( 19, J, K) 

V(3, J,K)=V(19, J,K) 

1/(3, J, SO =1/(19, 0,(0 
U(20, J,IC)=U(4, J,K) 

VC20,J,JO=V(4, J,IC) 

WC 20, 0,IO=:/C4, J,K) 

UC21, J,fC)=UC5, J, 10 
VC21, J,K)=V(5,J, K) 

W£ 21, J, l() =5/(5, J, K) 

U(22,.',K)=U(6, J,IC) 

V(22, J,K)=V(6, J,K) 

1/(22, J,!0=W(6, J,K) 

24 CONTINUE 

DO 22 1=1,22 

U( I , 1, X)=tJ( 1 , 17, K) 

V( I , 1, \ ) = V( I , 17, K) 

./ (I, 1, it )=!/((, 17, K) 

U( I , 2, 10 =U( I , 18, K) 

V( I , 2, K) =V( 1 , 18, fO 
V/C I, 2, K) =1/(|, 18, K) 

U( I , 3, 1C) =U( 1 , 19, K) 

VC I , 3, !C) =VC 1 , 19, 1C) 
i/( 1 , 3, K) =l/( | , 19, K) 

U( I , 20, IC)=U( 1 , 4, K) 

VC J , 20, IC)=V( I , 4, K) 

;/C 1 , 20 , 10 = 1 /( 1 , 4,10 

0(1,21,10=0(1,5,10 

VCl,21,IO=V(i,5,K) 

'/( I , 21, K) =W( 1 , 5, K) 

U(I,22,K)»U(I,6,K) 

V(l, 22,10=70,6,10 
’./(I, 22,10=1/(1,6,10 
22 CONTINUE 
20 CONTINUE 
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COMPUTE THE EDDY VISCOSITY COEFFICIENT 
no 25 K*2,6 
DO 25 J=2,21 
DO 25 ) “2, 21 
KK=K”1 
JJ=J-1 
1 1-1-1 

EV{ I I, JJ,KI0=C3*SQRTC2.*(U(1,J,K)-U(I-1,J,K))**2 
+2.*C V£ I , J, 10 -VC I , J-l, JO) **2+2.* CUC I, J,IC)-WC !, J,K-1))**2 

+ . 25*c cun , j+i, io-uc i , j, io+vc i + i, j, io-vc i , j, k) ) **2 

+ CUC I -1, J+1,K)~U( I -1, J,K)+V( I, J,K)-V(I-1, J,K))**2 
+ CUC I -1, J,K)-U(1-1, J-1,K)+V( 1 , J-1,K)-V( 1-1, J-l, K)}**2 
+ (UC l , J,K)-U( I , J-l, K)+V( I +1, J-l, K)-V( I , J-1,K) ) **2 
+ CVCI, J,K+1)-V(1, J,K)+W(I, J+1,K)-W(I, J,K))**2 
+ (7(1, J,K)-V(I,J, K-D+WC 1 , J+l, K-D-WC t , J, K-l) ) **2 
+ CVC 1, J-l, K) -VC I, J-l,K-l)+UCl , J, K-D-V/C I, J-l,K-l))+*2 
+ ( VC I , J-l, IC+D-VC I , J-l, 10 +WC 1 , J, to -W( 1 , J-l, 10 ) **2 
+ CW( I +1, J, 10 -WC I , J, lO+UC ! , J, K+l) -U C 1 , J, K) ) **2 
+CWCI+1, J,K-1 >-W( I, J,K-1)+U(I, J,lO-U( I, J,K-U)**2 
+ CWC I , J, K-D-WC l -1, J,K-D+UCl-l, J,K)-U( I -1, J,K-1))**2 
+(W( I , J,K)-WC I -1, J,K)+Ut 1 -1, J,K+1)-U( I -1, J- K) )**2 ) ) 

CONTINUE 

COMPUTE R.H.S. OF MOMENTUM EQUATION (PRESSURE GRADIENT 
NOT INCLUDED) 

IFUS.EQ.l) GO TO 26 
K=4 

uO TO 27 
DO 35 K=2, 4 
DO 35 J=2, 19 
DO 35 1=2,19 
KK=K-1 
JJ=J-1 
11*1-1 
Jl=l 
J2=l 

IF(dJ.NE.l) GO TO 58 

Jl=-1 

J2=-l 

S5X= 

-C4*C<U(I+2, J+1,K+D+U< |+1, J+1,K+1)>**2 
-CUCI+1, J+l, K+D+UC I , J+1,K+1))**2)-C5* 

(CUC 1+2, J+l, K+D+UC 1 +2, J+2, K+l)+U( 1+1, J+2, K+D+UC I +1, J+I,K+1>)**2 
-CU( I +1, J+l, K+D+UC I +1, J+2, K+l) +U( I , J+2, K+D+UC I, J+l, K+l) )**2 
+ (U( ! +2, J, K+l )+U( I +2, J+1,K+1)+U( 1+1, J+l, K+D+UC 1 + 1, J, K+l) )**2 
-CUC I +1, J, K+l)+U( 1+1, J + l, K+l ) +U ( I , J+l, K+l )+U( I , J,K+1))**2 
+ CUC 1+2, J+l,K+2)+UCI+2, J+1,K+1)+UC t+1, J+1,K+1)+U( 1+1, J+l,K+2})*+2 
-CUC I +1, J+l, K+2)+U( 1+1, J+1,K+1)+U( I , J+l, K+D+UC I , J+l, K+2 ) )**2 
+ (IJ{ I +2, J + l, K+D+UC 1+2, J+l, K)+U( 1+1, J+l, 10+UCI +1, J+l, K+l))** 2 
-CUC 1+1, J+l, K+D+UC l +1, J+l, IO+UC I , J+l, !0+U(T , J+l, K+l) )**2) 
SSX=$’SX-C4* 

' ( (UC t +1, J+1,K+1)+U ( 1 +1, J+2, K+l) )* (VC I +2, J+l, K+D+VC ! +1, J+l, K+l) ) 

-CUC I +1, J, K+D+UC 1+1, J+l, K+l) ) * C V C 1+2, J,1C+1) + V( I +1, J, K+l)) )-C6 
*C cue 1+2, J+l, K+D+UC! +2, J+2, K+D+UC 1+1, J+2, K+D+UC t+1, J+l, IC+U)* 
VC I +2, J+l, K+l) 

-CUC! +2, J, K+D+UC! +2, J+l, K+l )+U Cl +1, J+l, K+D+UC I +1, J, K+l))* 
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* VC1+2,J,K+1) 

* +(U( 1+1, J+1,K+1)+U( !+l, J+2,K+1)+UC I, J+2,K+D+U( l, J+l, K+D)* 

* VC! +1, J+l, K+1) 

* -cuu+i, j,k+d+uci+i, j+i,k+d+uci, j+i,k+d+uci, j,k+u>* 

* V(I+1,J,K+1))-C5* 

* (CUC 1+1, J+1,K+1)+UC 1+1, J+2,K+l)+U(!+l,J+l,K+2)+UCl+l, J+2, K+2))* 

* C VC I +2,J+1,K+1)+V(I+1,J+1,K+1)+VC [+2,J+l,K+2)+V( t +1, J+l, K+2 } ) 

* -CUC I +1, J, K+D+l/C J+l, J+1,K + D+U( I +1, J,K+2)+U( I +1, J+l,K+2) )* 

* (VC 1+2, J,K+D+VC I + I,d,K+l) + VC I +2, J,K+2)+V( 1 +1, J,K+2)) 

* +(UC 1+1, J+1J0+UC I +1, J+2,K)+UC 1+1, J+l, K+D+UC I +1, J+2, K+l))* 

* CV( I +2, J+1,1C) + V( 1 +1, J+1,K)+VC I +2, J+1,K+1)+V( [+1, J+l, K+D) 

* -CUC 1+1, J,K)+U( t +1, J+1,K)+WC 1 + 1, J, K+l)+U( 1+1, J+l, K+D )* 

* (VCl+2, J,lO + V( 1+1, J,IC)+V( !+2, J,K+1)+VC i+l,J,K+l)D 
SSX=SSX-C4 

* *( CUC 1+1, J+l, K+1J+UC I +1, d+l,K+2) )*(W< 1+2, d+1, K+D+WC 1 + 1, J+l, K+1)) 

* -CUC I +1, J+l, 10 +UC 1+1, d+1, K+l )>*(W( 1 +2, J+1,K )+!•/( t +1, d+1, K) ) }-C6 

* *((U(|+2, J+l, K+2) +U() +2, J + l, K+D + UC ! + l, J+l, K+D+UC I +1, J+l, K+2) ) * 

* l/C I +2, d+1, K+l) 

* -CUC I +2,d+l,K+l)+UC 1+2, J+1,K)+U( l+l,J+l,K)+UC1+l,J+l,K+l))* 

* WC I +2, d+1, K) 

* +CUC ! +1, d+1, K+2 )+U( 1 +1, J+l, K+D+UC I , J+l, K+l) +U C I , J+l,K+2))* 

* W( 1+1, d+1, K+l) 

* -CUC I +1, J+1,K+1)+UCI+1, J+l, K)+U( I , J+l, K)+U( 1, J+l, K+l))* 

* WCI+1, J+1,K))-C5* 

* ( CUC I +1, J+1,K+1)+U( t+1, J+l, K+2 ) +UC I +1, J+2,K+1)+UC 1+1, J+2,K+2))* 

* CWC I +2, J+l, K+D+l/C 1 +1, J+1,K+1)+WC f +2, J+2, K+l) +W< I +1, J+ 2, K+l) ) 

* -CUC I +1, J+l, K) +U( I +1, J+l, K+ 1 ) +U ( ! +1, J+2,K)+UC 1+1, J+2, K+l))* 

* CWC ! +2, J + l, K)+W( I +1, J+l/O+WC 1 +2, J*2,K)+t/C 1+1, J+2,K) ) 

* +CUC I +1, J, K+D+UC I +1, J,K+2)+UC 1 +1, J + 1,K+1)+UC I + 1, J+l, K + 2 ) )* 

* CWC I +.2, J, K+D+WC I +1, J, K+D+l/C I +2, ,1+1, K+l)*l/( I +1, d+1, K+l) } 

* -CUC I + 1, J, K)+U( ] +1, J, K+l )+UC 1 +1, d + 1, K)+U( t +1, J+l, K+D ) * 

* 1 + 2, J, K)+W( 1 + 1, J, K) +WC 1 + 2, J + l, K) +W (I+1,J+1,K))) 

5XC 1 1 , JJ,KK)=SSX 

* +(8.*(EVC 1+1, J,K)*(U([+2,d+l,K+l)-U( 1+1, J+1,K+1)) 

* -EVC I,d,K)*CU( 1+1, J+l,i;+l)-UCt , J+l, K+l))) 

* + CEVC1, J,K)+EVCI+1, J,K)+EVCI+1,J+1,K)+EVCI,J+1,K))* 

* CUC !+l,d+2,K+l)-UCl +1, J+l, K+l) +VC 1 + 2, J+l, K+D - VC f + 1, J+l, K+l) ) 

* - C EVC f ,d,K) + EV(I+l, J,K) + EV( I +1, d- 1, 10 + EVC I , d-l,K) )* 

* CUC I +1, J+1,K+1)-U( 1 +1, J,K+1)+V( 1+2, J, K+l) -VC 1 + 1 , J, K+1) ) 

* +CEVC !,d,K) + EVCl+l,d,K) + EVCl, J,K+1) + EVCI +1, J,K+D)* 

* CUC I +1, d+1, K+2)-UC I +1, d+1, K + D+WC I +2, J+l, K+D-WC I + 1, J+l, K+1)) 

* -CEVC1,J,K) + EVC!+1,J,K) + EVC1,J,K-D + EVCI+1,J,K-D)* 

* CUC I +1, J+l, K+U-UC I +1, J + l, K)+WC [+2, J+l, K) -W( 1 +1 , J+l, K) ) ) *C7 
GXC I I , JJ, KK)-SX C II , JJ, KK) 

* -2 . *Ud* ( C FLGATC JJ-1) *il- C . 5+11) ) *( (UC I + 2, J+l, K+1) - UC I , J+l, K+ U ) *C11 

* +CUC I + 2, J+l, K+2 ) + 2 . *U ( I +2 , J+l, K+D - 2, *UCI , J+l, K+1)-U(i , J+l, K+2) 

* +UC 1+2, J+1,K)-U( I , J+l, K) ) *C12 ) 

* +CFLQATC JJ-l)*ll+Q. 5*II-( .5+11) )* 

* ( U C I + 2, J+l, K+1) +UC I +2, J+2, K+D-UC 1 , J+2, K+1) -UC I , J+l, K+l) ) *C12 

* +CFLOATC Jd-D*H-0. 5*H-C .5 + 1)))* 

* CUC 1 +2, J, K+D +UC 1+2, J+l, K+D-UC t, J+l, K+D -UC 1 , J, K+l) )*C12 ) 

GXCI ! , JJ, KK) =SXC I J , JJ, KK) 

* -UO*CC CFL0ATCJJ-l)*‘ri+.5*!l-C . 5+U) ) * C VC I +2, J+l, IC+D+VC I + 1, J+ 1, K+l) ) 

* -CFLOATC JJ-1)*H-. 5*il-C .S+M))*(VC I +2, J, K+D+VC I +1, J,K+1) ) )*C11 

* +( CFLOATC JJ-D *11 + . 5*11- ( . 5+11) }*(i. *V( 1+2, J + l, K+l) 
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* -{FLOAT C JJ-1)*H-. 5*H-( . 5+11) ) *4 . *VC I +2, J, K+1) 

* + CFLOATC J J-l)*M+.5*tl-C . 5+li) )*4 . *V{ I +1, J+l, K+l) 

* -CFLOATC JJ-1) *11- . 5*H-( . 5+H } ) *4 . *V( I + l,d,K+l) 

* +CFLOATC JJ-1) *H+.5*H- ( . 5+H ) )* 

* (VC 1 + 2, J+l, K+l)+V( I +1, J+l, K+D+VC I + 2, J+l, K+2 )+ VC 1 +1, J+l, K+2) ) 

* -CFL0ATCJd-l)*H-.5*H-{ ,5+il) )* 

* C VC I +2/ J, K+D +V C !+l,J,K+l)+VCl+2, J, K+2)+V C 1 +1, J,K+2)> 

* + {FLOAT ( JJ-1) *11+ .5*11- { . 5+H ) )* 

* (V( 1+2, J+1,K) + V(I+1, J+l, K)+VU +2, J+l, K+D+VC 1 +1, J+l, K+1)) 

* -(FLOATC JJ-l)*ll-.5*H-( . 5+H) )* 

* (VC 1 + 2, J, IO+VC 1+1, J,K)+VC 1 +2, J, K+D + VC 1+1, J, K+D) )*C12) 

SXCt I , JJ, KK) =SX( I 1, JJ,KK) 

* -UQ*C(FLOAT( JJ-1) *H-( .5+11) )* 

* CUC I + 2, J+l, K+D+WC 1 +1, J+l, K+l )-W( 1 +2, J+l, K) - W( I + 1, J+l, K) )* 

* (C11+C2 ) 

* +C(FL0AT(JJ-l)*H+.5*H-£ .5+H))* 

* (WC 1 +2, J+l, K+D+WC I +1, J+l, K+l) +W( 1 + 2, J+2, K+D+UC J +1, J+2, K+l) 

* -W{ 1 +2,0 + 1, K) -W( 1 +1, 0 + 1, 10 - W{ t +2, 0+2, IO-WC t +1, J+2, K) ) 

* + CFL0ATC JJ-1 ) *H-. 5*H- (. 5+H) )* 

* OK I +2, 0, K+l) +l/{ 1 +1, 0, K+1)+VK I +2, J+l, K+1) +':/£ I + 1, J+l, K+l) 

* -WC I +2, 0,10 -VK 1 +1, 0, K)-W C 1 +2, 0+l,K)-U( I +1, 0+1, K) ) )*C12 ) 

* +2 . *0flEGA* . 25* 

* CV( 1+2, 0, K+l) +V{ ! +2, J+l, K+l) +VC I +1, 0+1, K+l)+V£ I +1, J, K+l) ) 

SSV--C4* 

* ({VC 1 +2,0+1, K+l )+V( I +1, 0+1, K+l) )*{UC I +1, 0 + l,K+l )+U{ I +1, J+2,K+1) ) 

* -CV( ! +1,0+1,K+1)+V( I, 0+1, K + 1)XU( I, J+l,K+D+l'U, J+2,K+1)))-C6 

* *C (V( I +2 , 0+1, K+D+VC 1+2, 0+2 , K+l) +V( ! +1, 0+2, K+l ) + VC t +1, 0+1, K+l) ) * 

* U{ ! +1, 0+2, K+l) 

* -(V( 1+1, 0+1, K+D+VC 1+1, 0+2, K+D + VC I, J+2, K+l)+V( I , J+l, K+l))* 

* UU, J+2, K+l) 

* +CVC1 +2,0, K+D+VC I + 2, J+l, K+l)+V( I +1,0+1, K+D+VC 1 + 1, J, K+l))* 

* UC1+1, J+l, K+l) 

* -CVC 1+1, J, K+D+VC ! +1,0+1, K+D+VC I, 0+1, K+D+VC I, J, K+l))* 

* UC 1 , J+l, K+l) )-C5* 

* C CVC I +2, J+l, K+D + V( I +1, J+l, K+l)+V( 1 +2, J + l, K+2J+VC 1 +1, J+l,K+2 ) )* 

* CUC 1+1, J+l,K+l)+U(l+l,d+2,K+l)+U([+l,J+l,K+2)+UCl+l,J+2,K+2)) 

* -CVC 1+1, J+l, K+1)+VC 1,0+1, K+1)+VC I +1, J+l, K+2)+VC I , J + l, K+2) )* 

* CUC1 , J + 1,K+1)+UCI, J+2,K+l)+UCt,J+l,!02)+UCI, J+2, K+2)) 

* +CVC 1+2, J+1,K)+VCI +1, 0+1,10 + V( 1+2, J+l, K+l)+V( I +1,0+1, K+l))* 

* (UC I +1,0+1,10+11(1+1, J+2,K)+U(l + 1, J+1,K+1)+UC 1+1, J+2, K+l) ) 

* -CVC 1+1, J + l, K)+V( 1 ,0+1,10 + VC 1 +1, J + l, K+D+VC 1 , 0+1, K+l))* 

* CUC 1, J+l,lO+U{|,Ji-2,IO+U(i , J+l, K+D+UC 1,0+2, K+D)) 

SSY=SSY-C4* 

* C ( VC J+l, J+2, K+D+VC 1+1, J+l, K+D) **2 

* -CVC 1+1, J+l, K+D+VC 1+1, J, K+D) **2)-C5* 

* C CVC I +2, J+l, K+D+VC I +2, J+2, K+D+VC 1 + 1, J + 2, K+D+VC I +1, J+l, K+l) )**2 

* -CVC 1+2, J, K+D+VC 1+2, J+l, K+D+VC 1+1, J+l, K+D+VC 1 + 1, J,K+1))**2 

* +CVC I +1,0+1, K+D+VC I +1, J+2, K+D+VC !, J+2, K+D+VC ! , J+l, K+1) )**2 

* -CVC 1+1, J, K+D + VC 1+1, J+l, K+D+VC 1,0+1, K+D+VC I , J, K+D) **2 

* +CVC I +1, J+l, K+1)+VC 1+1, J+2, K+D + VC I +l,d + 2, K+2) +VC 1+1, J+l, K+2) )**2 

* -CVC 1+1,0, K+D + VC 1 +1, J+l, K+D+VC 1+1, J+l, K+2 )+V( 1+1, J, K+2)) **2 

* +CVCI+1, J+l, lO+VC 1+1, J+2, IO+VC 1+1, J+2, K+D+VC 1 + 1, J+l, K+D) **2 

* -CVC 1+1, J, IO + VC l+l, J+l, K)+V Cl +1, J+l, K+D+VC 1+1, J, K+D )**2> 
SSY=SSY-C4* 

* (CVC I + 1,0+1, K+D+VC 1+1, J+l, K+2) )*CW( ! +1, J + l, K + D+WC l+l, J+2, K+D) 
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* -CVC I +1, J+l, ft) +V( I +1, J+1,K+1) ) *(WC I +1, J+1,K) +W( I +1, J+2, K) ) )-C5 

* *C(V(I+1, J+l, K+D+VC l +1, J+l, K+2) + V( I + 2, J+l, K+D+VC I +2, J+l, K+2) )* 

* 07(1+1, J+l, K+D+UC l +1, J+2,K+D+W<l+2,J+l,K+l)+W(I + 2, J+2, K+l)) 

* - C VC 1+1, J+1,K)+V( I +1, J+l, K+l)+V( 1+2, J+l, K)+VC 1+2, J+l, K+l))* 

* (UC I +1, J+l, lO+WU +1, J+2, K) +WC I +2, J+l, K)+U( I +2, J+2,K) ) 

* +( V( I , J+l, K+l) +V( ] , J+l, It+2) +V( I +1, J+l, K+D+VC i+l, J+l, K+2))* 

* CUC ! , J+l, K+D+UC l, J+2,K+l)+t7C [+1, J+1,K+1)+17CI+1,J+2,K+D) 

* -eve i,j+i,io+vu, j+i, k+d+vc i+i,j*i,k}+v( i+i,j+i,k+ id* 

* OK I , J+1,K)+V7( f , J+2, K) +17(1+1, J+l, K)-+W( l+l,J+2,K)))-C6* 

* ((VC I + 1, J+l, K+D+VC 1+1, J+2, K+D+VC 1+1, J+2,K+2)+V( f+1, J+l,K+2))* 

* UC I+l, J+2, K+l) 

* -*(V( ! +1, J+l, K)+V( I +1, J+2,K)+VCI+1, J+2, K+1)+VC I + l, J+l, K+l))* 

* 17(1+1, J+2,K) 

* +(V( I+l, J,K+1)+V( I+l, J+l, K+l) + V{ I+l, J+l,K+2)+V( I+l, J, K+2))* 

* i7( I+l, J+l, K+l) 

* -(VC 1+1, J,K) + VC I+l, J+1,K)+VC I+l, J+l, K+D+VC I+l, J,K+1))* 

* uci+i, j+i,io) 

SY( I I, JJ,KK)=SSY 

* +( ( EV( I, J,K)+EV(I , J+1,K)+EV( * , J+l, K+l )+EV( I , J,K+1) )* 

* CVC I+l, J+l,K+2)-V( I +1, J+l, K+l) +17(1 +1, J+2, K+l )-W( 1 + 1, J+l, K+l)) 

* -( EV( I , J, K) + EVC I , J, K-l ) + EV( I , J+l, K-D + EVC I , J+l, K) )* 

* CVC I+l, J+l, K+l)- VC I+l, J+l, K)+W{ I+l, J+2, K) -17(1+1, J+l, K)) 

* +3. * ( EV( I , J+l, JO* (VC I+l, J+2, K + l) -VC I+l, J+l, K+l) ) 

* -EV( J , J, !0 *( VC 1 +1, J+l, K+l )- V( I +1, J, K+l ) ) ) 

* + ( EV( 1 ,J,IO + EV( I + 1,J,I0+EVC I + l,J+l,it) + EV( ] , J+1,K)>* 

* CVC 1+2, J+l, K+l) -VC I+l, J+l, K+l) +IJC I+l, J+2, K+l) -UC I+l, J+l, K+l)) 

* -CEVC I , J, K) + EV( I , J+l, K) + EV( ! - 1, J+l, K) + EVC I -1, J,K))* 

* (VC I+l, J+l, K+l)- VC I , J+l, K+D+UC I , J+2, K+1)~U( I , J+l, K+l) ) )*C7 
SYCI I , JJ, KIO-5YC I I, JJ,KI0 

* -U0*( C FLOAK JJ-1 )*!!- Jl* . 5* (1 .+J2*H) ) * 

* ( CVC I +2, J+l, K+l) -VC I, J+l, K+l) )*C11 

* +C2 . *V( I +2, J + l, K+D+VC I + 2, J+l, K+2 )-2. *V( I , J+l, K+l) -VC I , J+l, K+2 ) 

* +VC 1+2, J+l, K) -VC I , J+l,IO)*C12) 

* +( FLO ATC JJ-1 )*H+ . 5*11- Jl*. 5*( 1 . + J2*!I) ) * 

* (VC 1 +2, J+l, K+D + VC I +2, J+2, K+l) -VC I, J+2, K+l)- VC I, J+l, K+l)) *C12 

* + C F LOATC JJ-1) *It- . 5*11- Jl* . 5*(1 . +J2*H) ) * 

* (VC I +2, J, K+l)+V( 1+2, J+l, K+l) -VC I, J+l, K+l )- VC I, J,K+1))*C12) 

* -2 . *OMEGA*0 . 25* 

* (UC I +1, J+l, K+l) +UC I + 1, J+2, K+l)+U( I , J+2, K+D+UC I , J+l, K+l) ) 

SSZ=-C4* 

* (07(I+2,J+l,K + l)+L7Cl+l,J + l,K+D)*CU(l+l,J+l,K+2)+UCl+l,J + l,K+D) 

* -07(1+1, J+1,K+1)+W( I , J+l,K+l))*(U(l,J+l,K+2)+U(l,J+l,K+l)))-C5 

* *C07C I +2, J+l, K+l) +WC I + l, J+l, K+l ) +17 C ! +2, J + 2, K+l )+WC I +1, J+2, K+l) )* 

* (U( I +1, J + l, K+2) +U( I +1, J+l, K+l)+U( I +1, J+2, K+2 )+U( I+l, J+2, K+l)) 

* - 07 ( I +1, J+l, K+l) +17 ( I , J+ 1, K+l) +W( I +1, J+2, K+D+l/C I , J+2, K+l ) ) * 

* CUC I, J+l, K+2) +U( I, J+l, K+D+UC I, J+2, K+2) +U( I, J+2, K+l)) 

* +0/C I +2, J, K+D+'./C I + 1, J, K+D +’/C I + 2, J+l, K+l )+’.!( f + 1, J+l, K+l))* 

* (U( I +1, J, K+2 )+U( I+l, J, K+D+UC I + l, J+l, K+2) +U( I+l, J+l, K+D) 

* -07C ] + 1,J, K+D+UC I, J,K+1)+W(1 +1, J+l, K + D +17(1, J + l, K+l))* 

* (UC I; J,K+2)+U(l , J, K+D+UC I, J+l, K+2 )+U ( I , J+l, K+l) ) ) -C6* 

* C C UC 1+2, J+l,K+2)+l/(l+2, J + 1,K+1)+W([+1,J+1,K+1)+17C1 + 1, J+l, K+2))* 

* UC I +1, J+l, K+2) 

* -0/C l +1, J+l, K+2) +17(1+1, J+l, K+D+UC I, J+l, K+D+UC I, J+l, K+2))* 

* U( I, J+l, K+2) 

* +CI7C 1+2, J+l, K+D+17C I+2-, J+l, K)+W( I +1, J+l, K) +U( I + 1, J+l, K+l) )* 


101 


Gf I p^ jPAGB ® 
OF POOR QUALITY 



* U( 1+1, J+l, K+D 

* -CWC 1+1, J+l, K+l)+W( I +1, J+l, K)+W( 1 , J+1,K) +W( 1 , J+l, K+1) )* 

* U( I , J+l, K+l)) 

SSZ=SSZ-C4* 

* CCUCl+l, J+1,K+1)+UCI+1, J+2,K+1))*CV(I+1, J+1,K+1)+V( 1+1, J+l, K+2)) 

* -UKl+1, J,K+1)+UCI+1, J+l,K+l))*(VCl+l, J,K+1)+V( 1 + 1, J,K+2)))-C5 

* *((W(I+1, J+1,K+1)+UC 1+1, J+2, K+D +W( 1+2, J+l, K+l )+U(l +2, J+2, K+l))* 

* (VC l+i, J+l-,K+i)TVCl+l,J+I,K+2)+VU+2,J+l,K+l)+V( I + 2 , J+l , K+ 2 J > 

* -(WU+1, J, K+D+WC l+l, J+l, K+D+WC 1+2, J, K+D+WC I +2, J+l, K+l) )* 

* C VC 1+1/0, K+l) +VO+1, J,K+2)+V( 1+2, J, K+l) +VC t +2, J,K+2)) 

* + CWC l , J+l, K+D+WC l , J+2, K+l ) +WC I +1, J+l, K+1) +WC I +1, J+2, K+l) )* 

* C VC 1/ J+l,K+D+VCl, J+l/K+2) + V( 1+1/ J + l/K+D+VC 1+1/ J+l,K+2)) 

* -CWC 1/ J,U+l)+Vin, J+l/ K+D+l/C 1 + 1/ J/ K+D+WC 1+1, J+l/K+l))* 

* C VC I / J/K+D+VC i / J, K+2) + V( l +1, J/K+l) +V( l + l/J/K+2)) )-C6* 

* (CWC1+1, J+l, K+U+lK l+l, J+2, K+D+WC 1+1, J+2, K+2) +UC l+l, J+l, K+2) >* 

* VCt+l/ J+l/K+2) 

* -CWU+1/ J/K+D+WC 1+1, J+l, K+D+WC 1+1/ J+l,K+2)+WCl+l,J,K+2))* 

* VC l + l/J/K+2) 

* +CWC1+1/ J + l, K)+WC 1 + 1, J+2, K)+WCt+l, J+2, K+D+WC 1+1, J+l, K+l))* 

* Vtl +1/ J+l, K+l) 

* -CUC l+l, J/K)+WC 1 + 1, J+l, K)+UC 1+1, J+l, K+D+WC 1+1, J, K+l))* 

* VC 1 + 1, J/K+l)) 

SSZ=SSZ-C4* 

* C CWC 1 +1/ J+l/ K+D+WC 1 +1/ J+l, K+2)) **2- 

* CWC 1 + 1, J+1,K)+WC! +1, J+1,K+1))**2)-C5* 

* C CWC l+l, J + l, K+D+WC 1 +1, J+l, K+ 2 ) +W C I +2, J+l, K+D+WC 1 +2, J+l, K+2) )**2 

* -CWC 1 +1, J+l, K)+WC l +1, J + l, K+D+WC 1 +2, J+l, 10+WC 1 +2, J + l, K+l) ) **2 

* +CWC 1, J+l, K+l)+W( 1 /J+l/K+2 )+WC I +1, J+l, K+D+WC 1+1, J+l, K+2)) **2 

* -CWC l , J+l/lO+WC I , J+l, K+D+WC 1 +1/ J+l, K)+WC 1 +1, J+l /K+l) )**2 

* + CWC 1 + 1, J+l, K+D+WC 1+1, J+l, K+2) +WC l+l, J+2, K+D+WC l+l, J+2, K+2)) **2 

* -CWC 1+1, J+l, IO+WC l+l, J+l, K+D+WC 1+1, J+2, K)+VK 1+1, J+2,K+1))**2 

* +CWC 1 + 1/ J/K+D+WC 1 +1/ J, ii+2 )+W ( 1+1, J+l, K+D+WC 1 +1, J+l/K+2) )**2 

* -CWC l+l, J/lO+WU + 1, J/K+l)+W( I +1, J+1,K)+W(1+1, J+1,IC+1))**2) 

SZCl l, JJ,KK)=SSZ 

* +CCEVCI, J/K+D+EVC 1+1, J/K + D+EVC l+l, J/K)+EVCI, J/K))* 

* (W( I +2, J+l, K+l) -WC l+l, J+l/K+D+UC 1+1, J+l, K + 2 )-U( l+l, J+l, K+l)) 

* -CEVC1 -1, J/K+D + EVC !, J/K+D +EVC 1, J,K) + EVU-1, J,K))* 

* (WC l+l, J+l, K+D-WC 1, J+l/K+D+UC I , J+l,K+2)-UC 1, J+l/K+D) 

* +CEVC1, J/K+D + EVC l, J/K) + EVC 1, J + 1,I0+EV( I , J+l, K+l))* 

* (WC l+l, J+2, K+D-WC I *1, J+l, K+l) +V( l+l, J+l, K+2) -V( l+l, J+l, K+l)) 

* -CEVCI, J-l/K+D + EVC 1, J-1,K) + EVCI/ J/IO + EVC I, J/K+l))* 

* CWC 1+1, J + l, K+D-WC I +1, J, K+D+VC 1 +1, J, K+2)-V( 1 +1, J, K+D ) 

* +3. *( EVC l , J/K+D*(W( I +1, J + DK+2) -WC I +1, J+l, K+D ) 

* -EVC l , J, K) *(W( I +1, J+l, K+D-WC l +1, J+l, K) ) ) } *C7 
SZCl 1 , J J, KK) -5ZC 1 j, JJ/KK) 

* — HO * C C F LOAT t J J-D *11- C * 5+M) ) *{ (WC I +2, J+l, K+D-WC I , J+l, K+l) )*C11 

* + CWC 1 +2, J+l, K+2) +2. *UC 1+2, J+l,K+l)-2. *WC I , J+l, K+l )-W( 1 , J+l, K+2) 

* +WC 1 +2, J+l, K)-U( 1 , J+l, K) )*C12) 

* +CFL0.1TC JJ-l)*H + 0 . 5*11- C . 5+H) )* 

* (WC 1 +2, J+l, K+D+WC 1+2, J+2, K+D-WC 1 , J+l, K+D-WC l , J+2,K+1) )*C12 

* +( F LOATC JJ-l)*M-0. 5*11- C .5+1!))* 

* (WCl+2, J, K+D+WC! +2, J+l, K+D-WC I , J, K+D-WC 1 , J+l, K+1))*C12) 

CO NT I HUE 

ADVANCE SX/SY/SZ IN TIME 
K»2 
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SMALLI NCQI (1,1),U1Q,1,M), 25D ) 

SMALL! NCQRC1,1),RX1C1,1,M), 25G) 

DO GO J=l,16 
DO GO E “1 ,16 

Ql (I, J)=QI Cl, J)+CALPHA*SXCl+l,J+l,K)-.5*QRCI, J))*DT 

60 CONTINUE 

SMALLOUT CQI ( 1, 1),U1( 1, 1,M), 256) 

DO 12 d*l, 16 
DO 12 1=1,16 
Ql (l,d)*SX(l+l, J+1,2) 

12 CONTINUE 

SMALLOUTCQI (1,1) ,RX1(1,1,M), 256) 

SMALL! NCQI (1, 1), VI ( 1, 1,M>, 256) 

SMALLI N(QR(1,1),RY1(1,1,M), 256) 

DO 61 J*l, 16 
DO 61 1=1,16 

Ql ( I, J)=QI Ct,d)+(ALPUA*SYC 1+1, d+l,K)-.5*QR( !, J))*DT 

61 CONTINUE 

SMALLOUTCQI Cl, 1), VI Cl, 1,M) , 256) 

DO 14 J-1,16 
DO 14 I =1,16 . 

Ql C I , J)=SYC I +1,0+1, 2) 

14 CONTINUE 

SMALLOUTCQI (1,1) ,RY1(1, 1,M), 256) 

SMALLI NCQI C1,1),W1(1, 1,M>, 256) 

SMALL ! MCQRC 1., 1) , RZ1C 1, 1,M) , 256) 

DO 62 J=l, 16 
DO 62 1=1,16 

Ql (I, J)=QI Cl, J) + CALPHA*SZ(I+1, J+l, I0-. 5*QR( I , J))*DT 
G2 CONTINUE 

SMALLOUTCQI Cl, 1),W1C1,1,M), 256) 

DO 15 d-1,16 
DO 15 1=1,16 
Ql C I , J)=SZ( I +1, J+l, 2) 

15 CONTINUE 

SMALLOUTCQI Cl, 1) ,RZ1C1,1,M), 256) 

C FORM THE R.rUS. OF POISSON EQUATION 

C1Q=1 ./ CDT*ALPHA) 

K=4 

DO 30 J=4, 19 
DO 30 1=4,19 
JJ=J-3 
11=1-3 

Ql C 1 1 , JJ) = CUC I , J, K)-U( I “1, J, IO + VC I , J, K)-VC 1 , 0-1, K)+!/C I , J, K) 

* -UC I, J,K-1))*C9 

QRC I 1, Jd)=QI Cl !,JJ)*C10+ 

* CSXC I -2, J-2, K-2 )-SXC I **3, J-2, K-2 ) +SYC I -2, J-2, K-2 )-3YC l-2,J-3,K-2)+ 

* S2C I -2, d-7, U-2J-SZC I ”2, d-2, IC-3 ) ) *C9 

C COMPUTE DIVERGENCE, TUROULENT ENERGY 

Dl V=D| V+QI Cl I, JJ) 

IJSQR=UC I , J, IO **2+USQR 
DUDX=(UCI, J,I0-UCI-1, J,K))Vll 
DUDX2=DUDX**2+DUDX2 
DUDX3 =DUDX**3+DUDX3 

EfJRG=ENRG+UC I , J, 10**2+ VC I , J, K)**2+W( I , J,K)**2 
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JO CONTINUE 

IF CM. EQ, 3) GO TO 9 
GO TO 7 

3 WRI TEC 6, 10J) 

103 FORMATC lil , 10X,*DI VERGENCE V *} 

WRITE (6/104) M 

104 FORMAT C 1H , 10X, ?.!l!!=, I 2) 

00 19 1=1,16/8 
WRITE (6,106) 1 

106 FORMATClil , lQX,2m=, 12) 

WRITEC6/105) (QJ(!,d},d=l,l6) 

105 FORMAT ( 2X, 3 ( E16. 7) ) 

19 CONTINUE 

7 SMALL0UT(QRU,1),PR1(1/1,M), 256) 

C SHIFT U, V,W,SX,SY,SZ ONE PLANE IN Z DIRECTION 

DO 32 K=1 , 6 
DO 32 J=l,22 
DO 32 1=1,22 
U( I , d/K) =U( I , J, IC+1) 

V(l,J,fO=V( J,d,K+l) 
l/(I,d,IO=W<l/d, K+l) 

32 CONTINUE 

DO 34 K=l, 2 

DO 34 d=l, 18 

DO 34 1=1,18 

SX{ I , J, K)“SX( I , d, l(+l) 

SY( I , d, K) =SY ( I , J,K+1) 

SZCI, J,tQ=SZ(l,J,K+l) 

34 CONTI NUE 

IFttl.EQ.lG) GO TO 38 
fl =11+1 
K2-K2+1 
GO TO 1 

38 DO 55 K=l, 16 
DO 55 d=l, 16 
DO 55 1=1/16 
PI 1 ( I / d/ K) =0. 0 
55 CONTINUE 

C INVERT THE POISSON EQUATION BY SUBROUTINE PSNEQN 

CALL PS ME ON 
WR I TEC 6, 110) 

110 FORMAT ( 1H ,10X,* PRESSURE FIELD *) 

SMALL I NCQRd/D/PRKl/l, 9)/ 256) 

' DO 47 1=1/16,8 

WRITECG/IOG) I 

WRITE(6,105) (QR(I/d)/d=l/16) 

47 CONTINUE 

C- ADVANCE PRESSURE GRADIENT IN TIME 

M=1 

K4=15 

50 DO 41 K=l, 3 

K3=K4+K 

IF(!(3.GT.16) K3=K3-16 

1 F(IC3 . GT. 16 ) IC3=K3-1B 

SMALL! N{QR(1,1), PR 1(1,1, K3), 256) 
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DO 42 0=1,16 
DO 42 1=1,16 
00= J+ 1 
11=1+1 

PR C I I , JJ,K)=QR( 1,0) 

42 CONTINUE 

C EXTEND PRESSURE FIELD PERIODICALLY IN X,Y DIRECTIONS 

DO 39 J=2,17 
PRC1, J,iO=PRC17, J,K) 

PRC 18, 0, K) = PR(2, J, K) 

39 CONTINUE 

DO 43 1=1,18 

PRC I , 1, K) =PR( E , 17, K) 

PRC 1 , 18, l() = PR( I , 2, IO 

43 CONTINUE 
41 CONTINUE 

K=2 

DO 54 J=l, 16 
DO 54 1*1,16 

QRC I , J)=-CPRC1 +2, 0+l,K)-PRCI+l, J+1,K))*C9 
54 CONTINUE 

SMALL! N(Q! (1, 1),U1C1,1,N ), 25G) 

DO 52 0=1, 1G 
DO 52 1=1,16 

Ql <1 ,J)=QI (I, J)+QP.(], J)*ALPiiA*DT 

52 CONTI NUE 

SMALLOUTCQI C 1, 1) , U1U, 1,M ), 256) 

SMALL t N(Qt ( 1, 1), RXK 1, 1, M ), 256) 

DO 53 0=1,16 
DO 53 1=1,16 
Qi Cl, J)=aiCl / O)+QRC!,0) 

53 CONTINUE 

SMALLOUTCQI Cl, 1),RX1C1,1,M ), 256) 

DO 64 0=1,16 
DO 64 1=1,16 

QRC i, J)==CPR Cl +1,0+2, K) -PRC 1 + 1, J+1,K))*C9 

64 CONTINUE 

SMALL I NCQt ( 1, 1) , VIC 1, 1,H ), 256) 

DO 65 1=1,16 
DO 65 0=1,16 

Qi (l,OJ=qi CI,0) + QII( I, J)*ALPUA*DT 

65 CONTINUE 

SMALLOUTCQI Cl, 1),V1C1,1,M ), 256) 

3MALL1HCQI C1,1),RY1C1,1,M ), 256) 

DO 63 0=1,16 
DO 6 3 1=1,16 
•HCI,0)=QIC1,J>+CLR<1,0) 

63 CONTINUE 

SMALLOUTCQI (1, 1),RY1C1,:.,M ), 256) 

DO 74 0 = 1,16 
DO 74 1=1,16 

QRC I ,0) = -( PRC t +1,0+1, K+D-PR Cl + 1,0+1, K))*C9 
74 CONTINUE 

SMALL I N(QI Cl, 1),W1C1, 1,M ), 25G) 

DO 75 1*1,16* 
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DO 75 J-1,16 

Ql (I, J)=QI ( I , J)+QR( I , J)*ALPHA*DT 
75 CONTINUE 

SMALLOJJTCQI C1,1),W1C1,1,M ), 256) 

SMALLINCQI C1,1),RZ1{1,1,M ), 256) 

DO 73 J«l, 16 
DO 73 1*1,16 
Ql { I , J)=QI Cl, d)+QR(I,J) 

73 CONTINUE 

SMALLOUTCQI U,1),RZ1U,1,M ), 256) 

I FCM, EQ. 16 ) GO TO 40 
M*M+1 
K4=K4+1 
GO TO 50 
40 L=L+1 

C COMPUTE TAYLOR MICROSCALE, SKEWNESS 

LAMBDA*SQRTC USQR/DUDX2 ) 

DUDX2=DUDX2/ 4096. 

DUDX3-DUDX3/ 4096. 

SKEVJ-UUDX3/DUDX2**!. 5 
ENRG*ENRG/ 4096. 

WR] TEC 6,111) SKEW, EMRG, D I V, LAMBDA 
111 FORMATC // , 10X, 1QHSKEWNESSX=, El 6. 7,5X, 7HENERGY*, E16.7, 
* 5X, Ilf 101 VERGENCE=, E16. 7, 5X, SULAML1DAX=, E1G. 7) 

D I V=0, 0 
ENRG=0. 0 
USQR=0.0 
DUDX2=0 . 0 
DUDX3=0.0 

IF(L.GT.O) ALPHA=1. 5 
IFCL.GT. 0) GO TO 78 
GO TO 6 
78 CONTI NUE 
STOP 
END 
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SUBROUTINE PSNEQM 

********** * ******************************************************** 


* POISSON SOLVER FOR 16 CUBE POINTS WITH PER I OH |C IJ.C. BY FFT * 

* ALGORITHM * 

* WAVE£I,J,K) IS THE WAVE NUMBER K' DEFINED IN CHAPTER III * 

* ISIGN=1 / I S I GN=-1 ARE FOR FORWARD AND BACKWARD TRANSFORMS, * 

* RESPECTIVELY * 


******************************************************************* 
LARGE PRlCifi, 16, 16 }, PI 1(16,16,16 ), WAVE (16, 16, 16) 

Dt MENS ION N£ 3 ), Ml ( G) ,M2 1 6), TRRC 16,3 5 , TRT (16, 3 ), E( 16) 

DIMENSION QR(16,16),QI (16/16) 

COMMON/SCM/QR/QI ,N,M1,M2,TRR,TRI ,E, I SIGN, PAI 
H=l./16. 

H2=H*H 
M3=H2*H 
NN-16 

NN3=NN*NN*NN 
C1-1./NN3 
C2-PAI /NN 
C3-1./H2 
I S ! GN“1 
CALL FFTX 
CALL FFTY 
CALL FFTZ. 

DO 26 K«l,16 
DO 26 J»l,16 
DO 26 1*1/16 

PR1U / J/K)=PR1CI / J,K)*C1/UAVE( 1/ J,K) 

PI l(t/J/K)~Pll(l/J / K)*C1/WAVE( 1 / 0/ K) 

26 CONTINUE 

PR1 ( 1/ 1/ 1) “0 . 0 
PI 1(1/1 / 1)=0.0 
I S 1 GN 3 ~1 
CALL FFTZ 
CALL FFTY 
CALL FFTX 
RETURN 
END 
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SUBROUTINE FFTX 

C ******************************************************************* 


c * fast Fourier transform in x direction * 

C * I SIGN®!/ I SI GN=-1 ARE FOR FORWARD AND BACKWARD TRANSFORMS, * 
C * RESPECTIVELY, THE FORWARD TRANSFORM RETURNS 16 TIMES * 
C * THE REQUIRED RESULT, FOR DEFINITION OF FPT SEE APPENDIX C * 


C ******************************************************************* 

LARGE PR1C16, 16,16), Pi I( IG, 16, 1C), WAVE ( 16,15,. 15} 

DIMENSION N( 3) ,M1C6) , M2 (6), TRR( 16, 3), TR 1(16, 3), E (16) 

DIMENSION QP.C1G, 16), Ql (16,16) 

COMMON/SCM/QR, Ql , N,M1,M2,TRR,TR1 , E, 1 S! GN, PAI 
DO 200 K“l, 16 

no 200 J=l,16 

DO 100 M«l, 3 
Jv)“N(M) 

KK-0 

5 I1=KK+JJ+1 
l(K“l l+JJ-1 
DO 10 I =1 I , KK 
L=1 -dJ 

Z=-PR1( L, J, K) 

T“-PI 1( L, J, K) 

PR1( L, J, K) — PR.1 £ I , J, K)+PR1( L, J,K) 

PI1(L, J,K)«Pll(l, J,K)+PI1(L, J,K) 

Z^PRIC 1 ,»J,K)+Z 
T=P1 1( I ,U,tO+T 

PR1(I , J,K)=Z*TRRU,M)-T*TRI (l,M)*ISIGN 
10 Pll( I, J,tO=Z*TRlCI ,M)*IS!GN+T*TRRCI ,M> 
l F (KK.NE. 16 ) GO TO 5 
100 CONTINUE 

DO 110 I a 2, 16,2 

PR1U,J,K)°-PR1(I,J,K)+PR1(1-1, J,K) 

110 PI 1( I ,U, K)=-PI 1(1 , J, K) + Pl 1( I ~1, J,K) 

DO 120 I “1, 15,2 

PR1( I , J, K) =*2, *PR1( I , J, K) -PR1( I +1, J,K) 

120 PI 1( I, J,K)*2,*PI 1(1, J, K)-P1 1(1 + 1, J, K) 

DO 20 L*l, 6 
I 1=M1( L) 

I 2=M2(L) 

DUM1-PR1(|2, J,K) 

PR1( I 2, J,K)”PRin 1, J,K) 

PRKIl, J,K)=DUM1 
DUM1=P! 1(12, J/K) 

PI 1( I 2,J,K)=PI 1( 1 1, J,K) 

Pll( II, J,I0=DUM1 
20 CONTINUE 
200 CONTINUE 
RETURN 
END' 
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SUBROUTINE FFTY 

C ******************************************************************* 

C * FAST FOURIER TRANSFORM IN Y DIRECTION * 

0 ******************************************************************* 

LARGE PR1C16,10,1G),P! If 1G, 16, 16 ) , WAVE{16, 16, 16 ) 

DIMENSION NC 35 ,111 C 6),M2 C 6 >, TRRf 16, 33 ,TRI C1G,3),E(16) 

DIMENSION QRC1G,1G),Q1 (16,16) 

COMMON/SCM/QR, Ql , N,M1, M2,TRR, TR1 , E, I S I GN, PA I 

DO 200 K=l, 16 

DO 200 1*1,16 

DO 100 M*l,3 

JJ=N(M) 

KK=Q 

5 I I “KK+JJ+l 
KK= t ! +JJ-1 
n o j-i Miff 
L=u-* J J 

Z=-PR1C I , L, K) 

T«-PI If J, L,K) 

PRlCI,L,fO=-PRl(l, J,K3+PRl(I,L,lO 
Pll(t,L,K)=PIlC!,.J,!C)+PH{f,L,lO 
Z=PRlf I, J,K)*Z 
T=P!1CI, J,K)*T 

PRK I , J, K)°Z*TRR( J,M}-T*TP. I (J,M)*ISIGN 
ID PI 1 ( I , J, K) =Z*TR I ( S! GN+T*TRR( «J,M) 

IFCKK.NE.16) GO TO 5 
100 CONTINUE 

DO 110 J*2,16,2 

PR1(I,J,K> — PRl(l,J,JO+PRl<|, 

110 Pi 1( I , J, lt)=-PI 1( I , <J, K) + PI 1( I , J-l, K) 

DO 120 0=1,15,2 

PR1C i , J,K) = 2i*PRl(I , J,IO-PRl(! , J+1,K) 

120 PI 1CI,J,K) = 2.*PI1(I,J, IO-PI1U, J+1,K) 

DO 20 L=l, 6 
!1=M1(U 
I 2=M2(L) 

DUM1=PR1CI, 12, K) 

PR1C I , 12, 10 =PR1C I , 1 1, K) 

PR 1C I , 1 1, K) =DUM1 
DUfll=P 1 1C I , I 2, K) 

PI1C1, 12,K)=PI1CI, I1,K) 

PI 1C I, I1,I0=DUM1 
20 CONTINUE 
20 Q CONTINUE 
RETURN 
END 
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SUBROUTINE FFTZ 

C ****v"t********************************ir***************************** 

C * FAST FOURIER TRANSFORM IN Z DIRECTION * 

C ***** ************************************************************** 

LARGE PR1U6,16,1G),PS 1 C16, 16 , 1G ) ,UAVE( 16 , 16, 16 ) 

DIMENSION NC3),M1(G),M2{6),TRR(1G,3),TR! (16,3),E<1G) 

DIMENSION QR(l6,16),ai (16,16) 

C0HM0N/SCM/QR,qi,M,Ml,M2,TRR,TRI # E,IS1GN,PAI 

DO 200 d*l,l6 

DO 200 1 «1, 16 

DO 100 M=l,3 

JJ=*N(M) 

KK*0 

5 1 1 -KK+JJ+1 

KK=l l+Jd-1 
DO 10 K=1!,KK 
L=K-JJ 

Z n -PR1( I , d, L) 

T=-P! 1( I , J, L) 

PR1C I , d, L) «PR1( I , d, K)+PR1( I , J, L) 

PI 1C I , d, L)=PI 1( I ,d,K) + PI K I , J, L) 

Z=PRK I ,d,K)+Z 
T*P I l(I,d,K)+T 

PRK I, d,K)=Z*TRR(K,M)-T*TRI CK,fU*tSlGN 
10 PI K I , d, K) -Z*TRI CfC,M)*fSIGN+T*TP,R(K,M> 

IFCKK.NE.16) GO TO 5 
100 CONTINUE 

DO 110 K a 2,16,2 

PRK l,d,IC)=-PRl(t , J,K)*PRKI, J,K-1) 

110 PlKi,d,K)=-PlKl,J,K) + PI l(!,d,K-l) 

DO 120 K a l, 15,2 

PRK I, J,I0=2,*PR1(I, J,K)-PRKI,d,K+l) 

120 PI 1( I, d,K)*2. *PIK!,J,K)-PI Kt ,d,K+13 
DO 20 L“l, 6 
I I«M1(L) 

I 2=M2 ( L) 

DUM1=PRKI,J, 12) 

PRKl,J,J2)=PRKl,J,ll) 

PRK l ,J, 1 1)=DUM1 
DUnl^P! KI,J, 12) 

PJ K I ,d, 1 2)-PI 1( I ,d, 1 1) 

PI K I ,d, I l)=DUMl 
20 CONTINUE 
200 CONTINUE 
RETURN 
END 
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SUBROUTINE TRRTR! 

************************************************************* *****£ 

* INITIATION of the trigonometric functions trr,tri required for * 

* THE FFT SUBROUTINES, AND WAVE (THE WAVE NUMBER K') FOR PSNEQN * 
******************************************************************* 

LARGE PR1( 16, 16,16), PI 1(16, 16, 16 ) ,WAVE{ 16,16, 16) 

DIMENSION N ( 3 ) , Ml (6),M2(6), TRR ( 16 , 3 ) , TR I (16,3),E(16) 

DIMENSION QRC 16,16), Ql (16,16) 

COMMOM/SCM/QR, Ql , N,Ml,M2, TRR,TRI , E, I SI GN, PA I 
DATA Ml/2,3,4,6,8,12/ 

DATA 1-12/9,5,13,11,15,14/ 

DATA E/0.,1.,2.,3,,4.,5.,6.,7.,-8.,-7.,-6.,-5.,-4.,-3.,-2.,-1./ 
11=1. / 1 6 * 

H2=*H*H 
NN“16 
C2=PAI/NN 
C3=l./H2 
M{ 1) “8 
N ( 2 ) “4 
N(3)=2 

DO 24 K*l, 16 
DO 24 0=1,16 
DO 24 1=1,16 

WAVEC1 , J,K)=-4.*(StN(C2*EC I ) )**2+SI M( C2*E( J) )**2 

* +SIN(C2*E(K))**2)*C3 
24 CONTINUE 

WAVE Cl, 1,1) =1.0 
DO 10 H=l,3 
L=N(M) 

J=1 

3Q K=J+2*M(M) -1 
DO 20 1 =J, K 

IFCI.GT.L) GO TO 22 
TPR( i ,M)=1 . 

TRI (l,W)«0. 

GO TO 20 

22 TRR(!,I1)=C0S(PAI*CI-1)/N(f1>) 

TRI ( I , M) =-S 1 N( PA I *( I "1 ) /N(M) ) 

20 CONTINUE 

IFCK.EQ.16) GO TO 10 
L=K+N(M) 

J=K+1 
GO TO 30 
10 CONTINUE 
RETURN 
END 
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SUBROUTINE INICON « 

C ******************************************************************* 

C * INITIATION OF A RANDOM FIELD WHICH SATISFY CONDITIONS SET IN * 

C * APPENDIX B. THE TWO REQUIRED FIELDS TO START WITH ARE * 

C * GENERATED IN SUBROUTINE RAMBM * 

C * SUBROUTINE SPECT3 EVALUATES THE THREE-DIMENSIONAL ENERGY * 


C * SPECTRUM OF THE GENERATED FIELD * 

C ******************************************************************* 

LARGE PR1(16,16,16),PI1(16,1G,16>,WAVEU6,16 # 16) 

LARGE U1C 16, 16, 19), VIC 16, 16, 19),W1(16,16, 19) 

LARGE U2(!6,16,18),V2(1G,16,18) / W2(16,16,18) 

DIMENSION QR(16,1G),QI £16^16) 

DIMENSION SR(16,16,6),SI t 10,16,6) 

DIMENSION N{ 3 ),M1( 6) , M2 ( 6 ) , TRR( 16 , 3 ), TRI (16,3),E(16) 

COMMON/SCM/QR, QI , N,M1,M2, TRR, TRI , E, I S I GN, PAI 
H e l. /16. 

HN-16. 

HN3=HN*HN*HN 
CALL RAN DM 
ARG»2.*PAI/HN 
DO 70 K«l, 9 
LL-K+9 

SMALL 1 N(SR(1,1,1),U1(1,1,K), 

SMALLINCSI £1,1,1), U1(1,1,LL 
SMALL I H(SR(1,1,2),V1£1,1, K), 

SMALL IN (SI (1,1, 2), VI Cl, 1, LL 
SMALL I tUSRC 1, 1, 3 ) , Wl(l, 1, 10 , 

SMALL I M( S 1 ( 1, 1, 3 ) , W1 C 1, 1, LL 
SMALL I N(SPv( 1, 1, 4), U2 £ 1, 1, K), 

SMALL! N(S I (1,1,4 ),U2£ 1,1, LL 
SMALL I N(SR( 1, 1, 5 ), V2( 1, 1> K) , 

SMALLI N(SI (1,1,5 ),V2 (1,1, LL 
SMALL! N(SR( 1,1,6 ),W2(1,1,K), 

SMALLINCSI ( 1,1,6 ),W2( 1,1, LL 
KK=*K-9 

DO 72 J*l,16 
DO 72 1-1,16 
JJ=J-1 
I 1*1 -1 

IFCJ.GT. 8) J J=J-17 
! F< I . GT . 8) 11=1-17 
ARG1=ARG* 1 1 
TRC1=1. -COS(ARGl) 

TRS1=SI N(ARGl) 

ARG l=ARG*dd 
TRC2=1.-C0S(ARG1) 

TRS2=S I N(ARGl) 

ARG1=ARG*KK 
TRC3=1 . “COS( ARG1) 

TRS3=Sf N(ARGl) 

A=TRC1*SR(I, d,4)+TRC2*SRCl, d, 5 ) +TRC3*SR( t , d,G) 

B*TRS1*SI (|,d,l)+TRS2*SI ( I , J, 2)+TRS3*SI ( l,d,3) . 

C=TRC1*S I C I , J, 4)+TRC2*Sl ( I ,.J, 5)+TRC3*S t (I , J,6) 

D=TRS 1*SR( I , J, 1 )+TRS2*SR( I , J, 2)+TRS3*SR( I , J, 3 ) 

E1=SR( t ,d,l)**2+SR( I ,d, 2)**2+SR( I , J,3 )**2 
E2=SRC»,d,4)**2+SP.( I , d, 5 ) **2+SR£ I , J,G)‘**2 


256) 

), 256) 
256 ) 

), 25G ) 
256) 

), 256) 
256) 

), 256) 
256) 

), 256) 
256) 

>, 256) 
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E3-SI (I , J,1)**2+SI C I , J, 2)**2+SI ( t ,J,3)**2 
E4=S I C 1 , J,4)**2+51 ( ! y J/ 5) **2+SI (!, Jy6)**2 

E5=SR( i / Jy 1)*SR( [ , Jy 4)+SRt I , J, 2)*SR( I , J, 5)+SR( i , J, 3)*SR( I , Jy 6) 
E6=S ! £ 1 ,0,13*51 ( I , J,4) + St ( I , J,2)*SI ( I, J,5)+SI C I, J,3)*S1 Cl, J,6) 
EEE=(El+E2+E3+E4)/2. 

DE)D“A*B*C*D 

t F (ODD . NE. 0 . 0) GO TO 71 
RD8 3 A*C 

l Ft B . EQ. 0.0. ANO . 8BB . NE. 0, 0) GO TO 75 
CCC*=B*D 

IFtC.EQ. 0.0. AND. CCC.NE. 0. 0) GO TO 77 
1 F( C, EG. 0.0. AND . D, EG, 0 . 0. AMD , B . HE, 0,0) GO TO 79 
l F C A. EQ. 0,0. AND . B . EQ. 0. 0 . ANO . C . ME. 0,0) GO TO 200 
! F( A, NE. 0. 0 . AND . B , EQ. 0, 0 , 4ND . C . EQ. 0.0. AND. 0. EQ, 0, 0 ) GO TO 210 
! F ( A. EQ. 0 . 0. AND. B . EG . 0. 0. AMD . C . EO,. 0. 0. AMD .0. ME. 0.0) GO TO 220 
IFtA.NE.0.0. AND . B. EQ. 0. 0, AND , C, EQ. 0. 0. AND .P.NE.0,0) GO TO 230 
I Ft A. EQ. 0 . 0 .AND. B. EQ. 0. 0. AND . C. EQ. 0. 0. AND. D. EQ. 0.0) GO TO 81 
71 RETA=RGEN(X1) 

i FCBETA .EQ. 0.0) GO TO 71 

771 GAMA=A/B*3ETA 
F1 5 =E1+D**2/C**2*E4 

F2=2. *BETA*ES-2. *GAf1A*D/C*EG 
F3=3£TA**2*E2+GAMA**2*E3~EEE 
F4=F2*F2-4. *F1*F3 
1FCF4.GT.O.O) GO TO 772 
BETA=BETA/2 , 

GO TO 771 

772 ALPHA=(-F2-SQRTCF4))/2./Fl 
DELTA=-D/C*ALPHA 

GO TO 73 

75 GAMA=P.GEM(X2) 

1 F( GAMA. EQ. 0.0) GO TO 75 
BETA“S/A*GAMA 
F1=E1+D**2/C**2*E4 
F2=2.*BETA*E5-2. *GAflA*D/C*E6 
F3 =3 ETA**2*E2+GAMA**2*E3-EEE 
F4=F2*F2-4,*F1*F3 
1FCF4.LT. 0.0) GO TO 75 
ALPHA=(-F2-SQRTCF4))/2./Fl 
DELTA=-D/C* ALPHA 
GO TO 73 

77 DELTA=RGENt X3 ) 

1 F (DELTA. EQ, 0.0) GO TO 77 
ALPHA="C/n*DELTA 
Fl-E2+A**2/S**2*E3 
F2=2. *ALPHA*E5 

F3=ALPHA**2*El+DELTA**2*E4+2 . *ALPHA*DELTA*E6-EEE 

F4=F2*F2-4. *F1*F3 

1FCF4.LT. 0.0) GO TO 77 

BETA 55 ( - F2-SQRTCF4) ) /2 . / FI 

GAI1A=BETA*A/B 

CO TO 73 

79 ALP1IA 55 RGENCX4) 

OELTA=nGEN(X5) 

IFCALPIIA.EQ. 0.0. OR. DELTA. EQ. 0.0) GO TO 79 
779 F1=E2*A**2/G**2*E3 
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F2=2.*ALPHA*E5 

F3=ALPMA**2*E1+0ELTA**2*E4+2,*ALPHA*DELTA*E6-EEE 

F4=F2*F2-4.*F1*F3 

1FCF4.GT. 0. 0) GO TO 780 

ALPHA«ALPHA/2. 

PELTA=DELTA/2. 

GO TO 773 

780 BETA=(-F2-SQRT(F4))/2./Fl 

GAMA=BETA*A/B 
GO TO 73 

200 BETA=RGEN(X6) 

GAMA=RGEN(X7) 

I FCBETA. EQ. 0. 0. OR, GAMA. EQ. 0. 0 ) GO TO 200 

201 F1=E1+D**2/C**2*E4 
F2=2.*CBETA*E5-GAMA*D/C*E6) 
F3=BETA**2*E2+GAMA**2*E3-EEE 
F4=F2*F2-4.*F1*F3 

I FCF4, GT. 0. 0) GO TO 202 
BETA-BETA/ 2. 

GAMA=GAMA/2. 

GO TO 201 

202 ALPMA=C-F2-SQRT(F4})/2./Fl 
DELTA=-P/C* ALPHA 

GO TO 73 

21G GAMA=RGENC X3 ) 

DELTA=RGEN(X9) 

BETA =0.0 

IFCGAMA.EQ. 0.0. OR. DELTA. EQ.O.O)GO TO 210 

211 F4 = C EEE“GAMA**2*E3“PELTA**2*E4-2. *GAMA*DELTA*EG )/El 
IFCF4.GT.0.0) GO TO 212 

GAI1A=GAMA/2. 

OELTA=DELTA/2. 

GO TO 211 

212 ALPHA=-SQRTCF4) 

GO TO 73 

220 GAMA=RGEM( X12) 

0ELTA=RGEN(X13) 

ALPHA=G . 0 

I FCGAHA . EQ. O.'O . OR. DELTA. EQ. Q . 0) GO TO 220 

221 F4 = ( EEE’*GAMA**2*E3-DELTA**2*E4-2. *GAMA*DELTA*E6 )/E2 
IPCF4.GT.0.0) GO TO 222 

G AMALGAM A/ 2. 

DELTA=DELTA/2. 

GO TO 221 

222 D£TA=-SQRT<F4) 

GO TO 73 

230 BETA=0 , 0 
ALP!1A“Q . 0 
GAHA=RGEM{ X14) 

I F ( GAMA. EQ. 0. 0) GO TO 230 
232 F1=E4 

F2=2, *GAMA*EG 
F3=GAMA**2*E3-EEE 
F4=F2*F2-4.*F1*F3 
1FCF4.GT, 0, 0) GO TO 231 
GAMA=GAMA/2. 



GO TO 232 

231 DELTA- < -F2-5QRT CF4))/2./Fl 
GO TO 73 
81 CONTINUE 
DELTA-0.0 
3ETA-Q.Q 
ALPHA-0. 0 
GAM A- 0.0 

73 SRC I , J, 1)=ALPHA*SP( I , J, 1) +BETA*SRC I , J, 4 } 
S I ( I , J, 1) =GAMA*S I C 1 , J, 1) +DELTA*S I C I , J, 4 ) 
SRC I , J, 2)=ALPHA*SRC I , «J, 2)+BETA*SRC 1,0,5) 
SI C I , J,2)=GAI1A*S1 C I , J, 2)+PE!.TA*S J C I , J,5) 
SRC I , J, 3)~ALPHA*SRC I , J, 3)+BETA*SRC I , J, 6) 
SI C I, J,3)=GAMA*S| Cl, J,3)+DE1.TA*SI C1,J,6) 

72 CONTINUE 

SMALLOUTC SRC 1,1,1), U1(1,1,K), 256) 
SMALLOUTCS! C 1, 1, 1 ) , U2C 1, 1, K) , 256) 
SMALLOUTC SRC 1, 1, 2 ), VI Cl, 1,K) , 256) 
SMALLOUTCS I C1,1,2),V2(1,1,K), 256) 
SMALL0UTCSRC1, 1, 3),W1C1,1, fO, 256) 
SMALLOUTCS! ( 1,1, 3 ),W2C 1,1,K), 256) 

70 CONTINUE 
ISIGN— 1 
DO 74 K=l,9 
M“K 

SMALL! N(QRC1,1),U1C1,1,K), 256) 

SMALL I NCQ! Cl, 1),U2(1, 1,K), 256) 

! FCM. EQ, 9) M-M-16 
MH«M+8 

SMALL0UTCQRC1,1),PR1C1,1,MM ), 256) 
SMALLOUTC Q I Cl,l),PIl(l, 1,MM ), 256) 

74 CONTINUE 

DO 76 K=10, 16 
LL-18-IC 
DO 76 J=l, 16 
DO 76 1=1,16 

IFC1 . EQ. 9.0R.J.EQ. 9) GO TO 1 
JJ=18- J 
I I =18-1 

I FC JJ. EQ, 17) JJ-1 
! F ( I 1 . E(l. 17 ) 11=1 
PR1C! 1 , JJ, LL) =PR1C I , J,K) 

PI It I I , J.l, LL) =-PllCI,J,K) 

GO TO 76 

1 PR 1C 1 , J, LL) =0, 0 
PI1C t,J,LL)=G.O 
76 CONTINUE 
CALL FFTZ 
CALL FFTY 
CALL FFTX 
00 78 K=l, 16 
no 78 J=l, 16 
DO 7S 1=1,16 
U1C I , J, K)=PR1C I ,*J,K) 

78 CONTINUE 

DO 84 K=l,9 
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M»K 

SHALL I N(QR(1, 1} , Vl(l, 1, K) , 256) 
SHALL! NCQI ( 1, 1 ) / V2C1, 1, K) , 256) 
IFCM.EQ, 9) M=M-16 
MM-M+8 

SMALL0UT(QRC1,1),PR1C1,1,MM ), 256) 
SMALLOUTCQI C1,1),PI 1(1,1, MM ), 256) 
8 A CONTINUE 

00 86 K=1Q,16 
LL-18-K 

no 86 J=l ,16 

DO 86 1=1,16 

I F( I .EQ. 9.0R.J.EQ, 9) GO TO 2 
J J— 18" J 
11*18-! 

I F (JJ . EQ. 17) JJ*1 
IFUl.Eq.17) ! I *1 
PR1( I I , JJ,LL)=PR1(1 , d,K) 

PI 1( ! I , JJ, LL)=-PI 1< I , J, K) 

GO TO 36 

2 PR1C I , J, LL)=0, 0 
PI 1CI , J,LL)=Q.O 

86 CONTINUE 
CALL FFTZ 
CALL FFTY 
CALL FFTX 
DO 88 K=l, 16 
DO 88 J=l, 16 
DO 83 1=1,16 
VKI, J,K)=PR1U,J,K) 

88 CONTINUE 

DO 94 K=l, 9 
M=K 

SMALL I fKQRC 1, 1) ,U1( 1, 1,K) , 256) 
SMALLIfKQI { 1, 1) ,W2( 1, 1, K) , 256) 
IFCM.EQ. 9) M=H"16 
HM=M+8 

SMALLOUTC QP.C 1, 1) , PR1C 1, 1, MM ), 256) 
SMALLOUT ( Q1 ( 1, 1) , PI 1( 1, 1, MM ), 256) 
94 CONTINUE 

DO 96 K“1G, 16 
LL=18-K 
DO 96 J=l, 16 
DO 96 1=1,16 

IFU.EQ. 9. OR. J. EQ. 9) GO TO 3 
JJ=13-J 

1 1=18-1 

IFfJJ.EQ. 17) JJ=1 
I F{ I t . EQ. 17) 11=1 
PR1CI I,JJ,LL)=PIU(I,J,K) 

PI 1( I I , JJ,LL)=-P! 1( l,d,K) 

GO TO 96 

3 PR1( I , J, LL) *0.0 
PI 1{ I , J, LL)=0. 0 

96 CONTINUE 
CALL FFTZ 
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CALL FFTY 
CALL FFTX 
DO 98 K-1,16 
DO 93 J«1,16 
DO 98 !«1,16 

van, j,k}*pri<i, j,k> 

98 CONTINUE 

CALL SPECT3 

RETURN 

END 


SUBROUTINE RAN DM 

******************************************************************* 
* INITIATION OF TWO RANDOM FIELDS REQUIRED FOR INICON * 

******************************************************************* 
LARGE PR1(1G,1G,16),PI1(16,1B,16),WAVE(16,16,16) 

LARGE U1 CIS, 10,19), V1C1G, 16, 19 ) ,W1C 16 , 16, 19 ) 

LARGE U2{ IS, 16, 18 ) , V2C 15, 16, 18 ), W2 C 16, 16, 18 ) 

DIMENSION QR(1S,16),QI (16,16) 

DIMENSION SR (16,16,0), SI (16,16,6) 

DIMENSION N(3),M1(6),M2(6),TRR(16,3),TR!(16,3),E(15) 

COMMON/SCM/QR, Ql , N,IIl,M2,TRR,TR I , E, IS! GN,PAI 
DIMENSION ENC14) 

-EN IS E(K), THE THREE-DIMENSIONAL ENERGY SPECTRUM, AND IS GIVEN 
-AS A FUNCTION OF K AT 2*PAI , 4*PAI, 6*PAI, ETC. 

DATA EN/. 52, 1.5,1. 3,. 96,. 66,. 46,. 33,. 23,. 16,. 105,. 072,. 05,. 032,. 0/ 
MN=16. 

ARG=2.*PAI/HN 
DO 20 K“l, 9 
KK=K-9 

DO 10 J«l, 16 
DO 10 ! “1, 16 
JJ=J-1 

IF( J.GT.8) JtJ=J-17 

I I «I -1 

I F( I . GT. 8 ) I 1*1-17 

I F ( I I . EQ. 0 .AND . JJ . EQ. 0. AMD. KK, EO„, 0) GO TO 10 
Y=RGEN( X) 

ZETA=2 . *PAI *Y 
Al=l.-COS(ARG*l I ) 

A2 = l. -COS( APvG*tlJ) 

A3=l . -COS( ARG*KK) 

R=SQRT(A1**2+A2**2+A3**2) 

PHI =ACOS ( A3/R) 

IFCA1.NE.0.) GO TO 5 
IFCA2.GT.O.O) THETA=PAt/2. 

IF(A2. LT.O.O) THETA=-PAI /2. 


oeiginalp^k 

POOR QUAINT* 


117 



5 


IFCA2.EQ. 0.0) THETA=2.*PAI*Y 
GO TO 2 
CONTINUE 

THETA=ATAN2CA2,A1) 

2 CONTINUE ' 

AMAG»SQRT(FLOAT( l 1*1 I +JJ*Jo+KK*KK) ) *2. *PAI 
WN=AHAG/ C 2 . *PA1 ) 

I VJf’M-m 
IWNP1=IUN+1 

ABSVEL= C ENC lt/N) + < EN< I WNP1) -ENC I WH) ) *CViN-IWH) )/C 2. *PAI *AMAG**2 ) 
Y1=RGEN(X1> 

VELR=SQRT( Y1*ABSVEL) 

VELI=SQRT(C1.-Y1)*ABSVEL) 

Ul( I , J,K> »-( COS(ZETA) *COS{ PHI )*COS(THETA)+S| N(ZETA)*S1 NCTHETA) ) 

* *VELR 

Vlt l , J,K) =( -COS (Z ETA)* COS C PHI 5 *S1 NCTHETA) +S I N(ZETA)*COS( THETA) ) 

* *VELR 

W1CI, J,K) =C0S CZETA) *8 I N( PH I ) *VELR 

UK I , J,K+9) =-{COSCZETA)*COS(PHI >*COS<THETA)+S I NCZETA)*St NCTHETA) > 

* *VELl 

VKl , J,K+9) = ( -COS(ZETA) *COS ( PHI )*SI N(THETA)+SI N(ZETA)*COS(THETA) ) 

* *VELl 

Win, J,K+9) =COS(ZETA)*S I N( PHI )*VELI 
10 CONTINUE 
20 CONTINUE 

Ul(l,l,9)“D.O 
Vl( 1, 1,9 )=0* 0 
U1C 1, 1, D ) =0. 0 
Ul( 1, 1, 18 ) =0. 0 
VIC 1, 1, 18) =0* 0 
UK 1, 1, 18 ) =0. 0 
00 GO K=l,9 
KK c K-9 

00 SO 0=1,16 
DO 50 1=1,16 
J J a J -*1 

IFCJ.GT.8) JJ=d-17 
11=1-1 

I F( I . GT . 8) I I =1 -27 

IF C I I . EQ.O.ANO.JJ.EQ. Q.ANP. KK. EQ. 0) GO JO 50 
Y=RGEN(X) 

ZETA=2.*PAl*Y 
A1=S | NC ARG* 1 1 ) 

A2=SIN(ARG*JJ) 

A3«*S1N{ARG*KK) 

R=SQRT( Al**2+A2**2+A3**2 ) 

PHI =ACOS(A3/R) 

IFCA1.NE.Q. ) GO TO 6 
1FCA2.GT.O.O) THETA=PA!/2, 

IFCA2.LT.0.U) TIIETA=-PAI/2. 

IFCA2. EQ.0.0) THETA=2.*PA1*Y 
GO TO 4 
6 CONTINUE 

THETA=ATAN2(A2, Al) 

4 CONTINUE 

At.iAG=SQRTC FLOAT (11 * I 1 +JJ*JJ+KK*KK) ) *2. *PA1 
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WN»AMAG/(2.*PAI ) 
fWN-WN 
I WNP1“! V/N+l 

ASSVEL*( EN( IWN)*(EH( I WNP1 ) -ENC I WN) ) * ( WN- 1 W» > ) / £ 2 , *PA 1 *AMAG**2 ) 
Y2=RGEN(X2) 

V£LR*SQRTCY2*ABSVEL) 

VELI=SaRT((l.-Y2)*ABSVEU 

U2CI, d,K) «“<C0S(2ETA>*C0S(PHI )*COSCTHETA)+SI N{ZETA)*SI NCTHETA)} 

* *VELR 

V2 U / J/ (0 “C-COS(ZETA)*COS{PHI )*S! NCTHETA) +S I N(ZETA) *C0S( THETA) ) 

* *VELR 

W2CI , J,K) =CaSCZETA)*SlMCPIH )*VELR 

IJ2( I , J, K+9 ) — CCOS(Z£TA)*COS(PHI )*C0S( THETA) +S I N(ZETA) *S I NCTHETA ) ) 

* *VELI 

V2 C I , J,K+9) = (-COS(ZETA)*COS( P!H ) *S t MCTHETA) +S I N(Z ETA) *C0S (THETA) ) 

* *VELI 

W2Cl,J,K+9) =COS(ZETA)*SI N( PHI )*VEL! 

50 CONTINUE 
60 CONTINUE 

U2 C I, 1, 9) *0. 0 

V2(l,l / n)«0.Q 

W2C1,1,9)*0.0 

U2 ( 1, 1, IS )«0, 0 

V2( 1, 1,18)=0.0 

W2C1,1,18)«0.0 

RETURN 

END 
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SUBROUTINE SPECT3 

C *************************************************** **************** 

C * EVALUATION OF THREE-DIMENSIONAL ENERGY SPECTRUM * 

C ******************************************************************* 

LARGE PR1(1G,16,1G),P11C1G,16,16), ENERGY (16,16,16} 

LARGE UK 16, 16, 19 ), VIC 1G, 16, 19), Wl{ 16, 16, 19) 

DIMENSION QRC16,1G),Q! <16,16) 

DIMENSION N( 3) ,M1(6) ,M2(6) ,TRR{ 16, 3) , TRI ( 16, 3), EC 16 ) 

COMMON/SCM/QR, C1I , N, Ml, M2, TRR, TR I , E, I S I GN, PAI 

PAI 2=PA1 *PAI 

H»l./16. 

H2=li*H 

H3=H2*H 

NN=16 

NN3-NN*NN*NN 
C=0. 6 

C TRANSFORMATION OF VELOCITY FIELD TO WAVE NUMBER SPACE, AND 

C FORMATION OF THE TURBULENT KINETIC ENERGY 

I S I GN=1 
DO 10 J»l, 16 
DO 10 1=1,16 
qiu,j)=o.o 
10 CONTINUE 

DO 12 IC=1, 16 

SMALL I NCQR(1,1),U1(1,1,K), 256) 

SMALLOUTC QR( 1, 1) , PRK 1, 1, K), 256) 

SMALLOUTCQI C 1,1), PI 1C1,1,K), 256) 

12 CONTINUE 
CALL FFTX 
CALL FFTY 
CALL FFTZ 
DO 14 K“l, 16 
DO 14 J=l, 16 
DO 14 1=1,16 

PRlC I , J, K)=CPR1C ( , J,fO /NN3)**2+C PI 1C 1 , J, KJ/NN3 )**2 
14 CONTINUE 

DO 16 K=l, 16 

SMALLI N CQRC 1, 1), PR1C1, 1, K), 256) 

SMALLOUTCQRC 1, 1), ENERGYC 1, 1,K), 256) 

16 CONTINUE 

DO 20 J=l, 16 
DO 20 1=1,16 
Ql Cl, J)=0.0 
20 CONTINUE 

DO 22 K**l,16 

SMALL! NCORC 1, 1), VIC 1, 1, K) , 256) 

SMALLOUT (QR( 1,1) , PR1 ( 1, 1, K) , 256) 

SMALLOUTCQI Cl, 1), PI 1(1,1, K), 256) 

22 CONTINUE 
CALL FFTX 
CALL FFTY 
CALL FFTZ 
DO 24 K=l, 16 
no 24 J«l,16 
DO 24 1=1,16 

PRlC I , J, K) =C PP.1C I , J, K) /NM3)**2 + ( PI 1C ! , J, IO/NN3 )**2 
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24 CONTINUE 

DO 26 K«l,16 

SMALLI NCQRC 1,1), PR1C1,1,K), 256) 

SMALL! NCQI (1, 1), ENERGY ( 1, 1,10, 256) 

DO 28 J»l,16 
DO 28 I 16 
QI C I,U3“Qf Cl , JJ+QRC I, J> 

28 CONTINUE 

SMALLOUTCQI C 1, 1) , ENERGY C 1, 1,K) , 256) 

26 CONTINUE 

DO 30 0*1,15 
DO 30 1=1,16 
QICI,J)«0.0 
30 CONTINUE 

DO 32 K-1,16 

SI1ALL1 NCQRC 1, 1) ,W1C 1, 1, K) , 256) 

SMALLOUTCQRC 1,1),PR1(1,1,K), 256) 

SMALLOUTCQI C 1,1), PI 1(1, 1,K), 256) 

32 CONTINUE 
CALL FFTX 
CALL FFTY 
CALL FFTZ 
DO 34 K-1,16 
DO 34 J=l, 16 
DO 34 1=1,16 

PR1CI, J,IC)-CPR1CI, J, K)/NN3)**2+CPI 1( I , J, K)/NN3)**2 
34 CONTINUE 

DO 36 K-1,16 

SMALL I NCQRC 1, 1),PR1C1, 1,K), 256} 

SMALL l NCQI Cl, 1), ENERGY C1,1,K), 256) 

DO 38 J-1,16 
DO 38 1-1,16 

QIC I,J)-(QI(I, J)+QRCI ,J)> 

38 CONTINUE 

SMALLOUTCQI (1,1), ENERGY ( 1, 1, K), 256) 

36 CONTINUE 

C EVALUATION OF THE 3-0 ENERGY SPECTRUM 

V/R I T E C 6 , 1 0 0 ) 

100 FORMAT C 111 , 10X, *TltREE-D I MENS I ONAL ENERGY SPECTRUM*) 
M-l 

EE-0 , 0 

WH-F LOATCM-l )*2. *PA! 

LL=1 

WR ITEC 5, 101) V/N, EE 

101 FORMATC 1H , 10X, 2IIK-, F6. 2, 5X, 5I!E{ K) =, E16. 7) 

DO 42 M-2,15 

WN-FLOAT CM-1)*2, *PAs 
EE-O.O , 

RAOP-WN+PA! 

RADM=UH-*PAI 

LL-0 

DO 44 K-1,16 

SMALL I NCQRC 1, 1) , ENERGYC 1, 1, K>, 256) 

DO 44 J-l, 16 
DO 44 1-1,16 
KK-K-1 
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IFCK.GT. 8) KK-K-17 

IFCJ.GT. 8) JJ-J-17 
I 1-1-1 

I F C I .GT. 8) 11-1-17 
R2=FL0AT(I 1*1 I +JJ*JJ + KK*KK> 

RR-SQRTt R2)*Z. *PAl 

I F ( RR. LT . RAD P . AND . RR * GE . RADM ) GO TO 45 
GO TO 45 
46 LL-LL+1 

EE-EE+QRC l / J) 

45 CONTINUE 
44 CONTINUE 

EE=EE/(FLOAT( LL) ) 

EE=4.*PAI*C*EE*WN*WN 
WRITEC >,101) WN, EE 
42 C0NT11UE 

CALL * RRTRI 

RETURI 

END 
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OF time step t* r> 

U CC"<*0*FNr nF vet *mr» 

I* i ... 

-?,«576*i 3E-02 *7.03u7339E-o3 2.l57i353E*02 S .«»? 01 19E-02 
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